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ABSTRACT 

Context. A first characterization of extrasolar planets by the observational determination of the radius has recently been achieved for 
a large number of planets. For some planets, a measurement of the luminosity has also been possible, with many more directly im- 
aged planets expected in the near future. The statistical characterization of exoplanets through their mass-radius and mass-luminosity 
diagram is becoming possible. This is for planet formation and evolution theory of similar importance as the mass-distance diagram. 
Aims. Our aim is to extend our planet formation model into a coupled formation and evolution model. We want to calculate from 
one single model in a self-consistent way all basic quantities describing a planet: its mass, semimajor axis, composition, radius and 
luminosity. We then want to use this model for population synthesis calculations. 

Methods. In this and a companion paper, we show how we solve the structure equations describing the gaseous envelope of a pro- 
toplanet not only during the early formation phase, but also during the gas runaway accretion phase, and during the evolutionary 
phase at constant mass on Gyr timescales. We improve the model further with a new prescription for the disk-limited gas accretion 
rate, an internal structure model for the planetary core assuming a differentiated interior, and the inclusion of radioactive decay as an 
additional heat source in the core. 

Results. We study the in situ formation and evolution of Jupiter, the mass-radius relationship of giant planets, the influence of the core 
mass on the radius and the luminosity both in the "hot start" and the "cold start" scenario. We put special emphasis on the validation 
of the model by comparison with other models of planet formation and evolution. We find that our results agree very well with those 
of more complex models, despite a number of simplifications we make in our calculations. 

Conclusions. The upgraded model yields the most important physical quantities describing a planet from its beginning as a tiny seed 
embryo to a Gyr old planet. This is the case for all planets in a synthetic planetary population. Therefore, we can now use self- 
consistently the observational constraints coming from all major observational techniques. This is important in a time where different 
techniques yield constraints on very diverse sub-populations of planets, and where its is difficult to put all these constraints together 
in one coherent picture. Our comprehensive formation and evolution model should be helpful in this situation for the understanding 
of exoplanets. 

Key words, stars: planetary systems - planet-disk interactions - planets and satellites: formation - planets and satellites: interiors - 
planets and satellites: individual: Jupiter - methods: numerical 



1. Introduction 

The number of known transiting extrasolar planets or planet 
candidates has recently increased exponentially, thanks both 
to ground-based observations (e.g. Gillon et al. 2011), and to 
space missions like CoRoT (e.g. Leger et al. 2009) and Kepler 
(Borucki et al. 2011). Combined with radial velocity measure- 
ments which yield the mass of the planet, one gets the planetary 
mass-radius (M-R) diagram, which is an observational result of 
similar importance as the semimajor axis-mass (a-M) diagram. 

The latter relation is available through the success of ongoing 
radial velocity surveys (e.g. Mayor et al. 2011). It is a goal of 
population synthesis models to understand the structure of the a- 
M distribution, due to the multitude of clues it contains for planet 
formation theory (e.g. Ida & Lin 2004; Mordasini et al. 2009a). 
A recent comparison of numerous theoretical and observational 
results mostly obtained by the radial velocity technique can be 
found in Alibert et al. (201 1) and Mordasini et al. (2012a). 
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The reason for the important role of the M-R diagram which 
is now available for a statistically significant number of planets 
is that one can derive the mean density of the planet. This con- 
strains the internal planetary structure which is of central impor- 
tance to understand the nature (Leconte et al. 201 1), but also the 
formation of the planet. The formation and evolution of the plan- 
etary mass-radius relationship is studied in the companion paper 
Mordasini et al. (2012b), hereafter Paper II (see also Mordasini 
etal. 2011c). 

Besides transiting planets, also the number of planets de- 
tected by direct imaging has increased significantly in the past 
few years, even though in absolute numbers, much less such 
planets have been found to date. But already these discoveries, 
like the planetary system around HR 8799 (Marois et al. 2008, 
2010) have triggered numerous theoretical studies regarding the 
formation of these objects (e.g. Dodson-Robinson et al. 2009; 
Rratter et al. 2010). Two points about these planets are particu- 
larly interesting: Their large semimajor axis a and the fact that 
we measure the luminosity L of young giant planets at some 
time t. Both quantities are important to understand the forma- 
tion mechanism (e.g. Marley et al. 2007; Janson et al. 2011). In 
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the near future, more capable instruments like SPHERE at the 
VLT (Beuzit et al. 2007) or GPI at Gemini South (McBride et 
al. 201 1) and later EPICS at the E-ELT (Kasper et al. 2008) wiU 
come online. We can therefore expect that the number of points 
we can put in the t-L, the a-L and (for cases with an independent, 
dynamical mass determination) the M-L diagram will increase 
from a handful at the moment to hundreds in a few years, sim- 
ilar to what has happened for the M-R diagram in the past few 
years. 

This shift from an era of discovery to one of a first physi- 
cal characterization of extrasolar planets by their radii and lu- 
minosities has profound implications for planet formation the- 
ory. Until recently, planet formation models studied mostly giant 
planets detected by radial velocity measurement, of which only 
the minimal mass and some orbital elements were known. Now 
we are confronted with a multitude of different sets of observa- 
tional data and constraints, each regarding primarily planets of 
different types: Transiting, close-in planets, most of them with 
a small radius (as found by Kepler, see also Paper II); directly 
imagined, self-luminous massive young giant planets at large 
distances from their host star; and a wide range in masses go- 
ing from super-Earth to Jovian mass planets at distances varying 
from close to the star out to several AU found by RV measure- 
ments (e.g. Mayor et al. 201 1). 

The final goal of planet formation theory to develop one the- 
ory to explain the formation (and evolution) of all these very 
different planets is a challenging one, and will require many ef- 
forts in the coming years. Nevertheless, the approach to bring 
together all these different observational data in a coherent way 
(which in itself is a non-trivial task, cf. Wolfgang & Laughlin 
2012) to constrain planet formation and evolution theories seems 
to be a promising route. In the end, one is not interested in a the- 
ory which can explain certain types of planets, but fails for other 
classes. 

In this and the companion paper we take a first step in this 
direction. We present multiple upgrades of our formation model 
(introduced first in Alibert, Mordasini & Benz 2004). The most 
important addition is that we now calculate not only the forma- 
tion of the planets, but couple formation in a self-consistent way 
with the subsequent evolution at constant mass once the proto- 
planetary disk is gone. With this approach, we can now calculate 
directly from one single model not only the planet's mass and 
semimajor axis, but also the planet's main physical characteris- 
tics like the radius, luminosity, surface gravity, effective temper- 
ature as well as the composition in terms of iron and silicates, 
ices and H.2/He. 

These basic characteristics are available for a planet at any 
time during its "life" starting as a tiny sub-Earth mass seed in 
the protoplanetary disk to a mature, billion of years old planet. 
A direct coupling of the planet's formation and its evolution is 
necessary, as it is well known that the formation has important, 
direct consequences for the evolution, in particular for the lu- 
minosity of giant planets at young ages ("cold" vs "hot" start 
models, cf. Marley et al. 2007; Spiegel & Burrows 2012). With 
this model development, we can now compare our simulations 
directly with observational constraints coming from radial ve- 
locity (and microlensing), transits as well as direct imaging. 

Other upgrades are the following: a detailed description for 
the rate at which gas is accreted by the planet in the disk-limited 
gas runaway accretion phase, an internal structure model for the 
solid (iron/silicate and possibly ice) part of the planet, the in- 
clusion of radioactive decay for the luminosity of the core, a 
new initial profile for the gaseous disk, a new prescription for 
the photoevaporation of the disk, including external and internal 



photoevaporation, and finally a realistically low grain opacity for 
the gaseous envelope (presented in a dedicated work, Mordasini 
etal. 2012c). 

We thus deal in this paper and in Paper II mostly with im- 
provements of the physical description of one planet. In other pa- 
pers, we addressed upgrades regarding the disk model (Fouchet 
et al. 2011), the migration of low-mass planets (Dittkrist et al. 
in prep) or the effect of the concurrent formation of several em- 
bryos in one disk (Alibert et al. in prep.). 

1.1. Organization of the paper 

The organization of the paper is as follows: In Sect. 2 we give a 
short overview of the model. In Section 3, we describe the modi- 
fications of the computational module that describes the gaseous 
envelope structure of the planet, extending it to calculate the 
structure not only during the pre-runaway formation phase as 
in our previous models, but also during the gas runaway accre- 
tion/collapse phase and the subsequent evolutionary phase af- 
ter the disk is gone. In Section 4 we address a related subject, 
namely how to calculate the gas accretion rate in the disk-limited 
regime, i.e. once gas runaway accretion of forming giant planets 
has started. 

In Paper II, we present further upgrades regarding the planet 
module, namely a realistic model for the density of the solid 
core of the planets and the inclusion of radiogenic heating in it. 
In Paper II we also describe shortly some modifications regard- 
ing the protoplanetary disk model. Finally, we use in Paper II 
the upgraded model in population synthesis calculations to study 
the formation and evolution of the planetary M-R diagram, the 
distribution of planetary radii, and the comparison with observa- 
tional data. 

In the remainder of this first paper, we show specific results 
obtained with the upgraded model: In Sect. 6 we study the cou- 
pled formation and evolution of a Jovian mass planet at 5.2 AU 
from the Sun. Many of the effects seen during this particular 
simulation are characteristic for the effects encountered during 
the formation and evolution of the planets in a general popula- 
tion synthesis calculation (Paper II). In Sect. 7 and 8 we discuss 
our results concerning the radii and luminosities of giant plan- 
ets, putting special weight on the comparison with other, more 
complex models. The luminosity of young Jupiters is further ad- 
dressed in a dedicated paper (Mordasini et al. in prep.). 

2. Combined model of planet formation and 
evolution 

The formation model used here relies on the core accretion 
paradigm, coupled to standard models of protoplanetary disk 
evolution and tidal migration of protoplanets. 

The basic concept of core accretion is to follow the concur- 
rent growth of an initially small solid core and its surrounding 
gaseous envelope, embedded in a disk of planetesimals and gas 
(Perri & Cameron 1974; Mizuno et al. 1978; Bodenheimer & 
Pollack 1986). Within the core accretion paradigm, giant planet 
formation happens as a two step process: first a solid core with 
a critical mass (of order 10 M @ ) must form, then the rapid ac- 
cretion of a massive gaseous envelope sets in. This process is 
thought to take typically several million years. Two other pro- 
cesses in the protoplanetary disk happen on a similar timescale: 
the evolution of the protoplanetary disk itself, and the orbital mi- 
gration of the protoplanets due to angular momentum exchange 
with the gaseous disk. As described in Alibert et al. (2005), we 
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have therefore coupled a classical core accretion model (very 
similar to Pollack et al. 1996) to models describing the latter two 
processes, using for the disk a standard a model, and prescrip- 
tions for the type I and II migration of the protoplanets. While we 
use simple, ID models for the individual processes, we consider 
the full coupling between them, which alone leads to complex 
formation scenarios (cf. Mordasini et al. 2009a). 

2.1. Limitations: no accretion after disk dissipation, 
primordial H 2 /He envelopes only 

From its conception, our model mainly deals with the formation 
and evolution of relatively massive, Neptunian and Jovian plan- 
ets even though that the vast majority of planets that actually 
form in a population synthesis are low-mass, failed core (proto- 
terrestrial) planets (Mordasini et al. 2009b). We currently do not 
include the giant impact stage (after the gaseous disk has disap- 
peared) during which terrestrial planet acquire their final mass. 
Therefore, our model is incomplete in the description of the for- 
mation of planets less massive than ~ 10M e (see Mordasini et 
al. 2009a for a discussion). 

Regarding the long-term evolution, we only consider primor- 
dial H_2/He envelopes for which we assume that the mass remains 
constant after the protoplanetary disk is gone, neglecting a pos- 
sible evaporation of the primordial gaseous envelope, as well 
as the outgassing of a secondary atmosphere. Such atmospheric 
mass loss for planets close to the host star will be considered 
in future work. Note that the minimal allowed semimajor axis 
for model planets is at the moment 0. 1 AU, so that no simulated 
planets are exposed to the very intense radiation field occurring 
on very tight orbits like for example CoRoT-7b (a = 0.017 AU), 
where atmospheric mass loss may be very important (Valencia 
et al. 2010). 

An exact lower boundary in mass down to which the model 
can be applied is difficult to specify, as we find in agreement 
with Rogers et al. (201 1) that already planets of just a few Earth 
masses can accrete non-negligible amounts of FL/He during the 
presence of the nebula, provided that the grain opacity in the 
envelope is low, which is probably the case (Movshovitz et al. 
2010). In Paper II we study the formation and evolution of a 
close-in super-Earth (~ 4M ffi ) planet with a tenuous 1% primor- 
dial H-2/He envelope. We find good agreement of the radius with 
the study of Rogers et al. (2011), showing that the evolution also 
of this kind of low -mass planets can be modeled. 

We describe in this paper only improvements (or the inclu- 
sion of new physical effects) relative to the original model de- 
scribed in Alibert et al. (2005), to which the reader is referred to 
for the remaining description. 



3. Gaseous envelope calculation during the 
attached, detached and evolutionary phase 

We now turn to the most important improvement of the model, 
which is the ability to calculate planetary envelope structures 
during the entire formation and subsequent evolution of the plan- 
ets. In previous models, we calculated the gaseous envelope only 
during the phase when the planet had a subcritical core mass less 
than the one needed to trigger gas runaway accretion (about 10 
M m ). This is sufficient if one is interested in the mass of the plan- 
ets only. In order to characterize the planetary structure (and thus 
to have R and L, too), we now include also the gas runaway phase 
(which is attained by planets becoming supercritical during the 
lifetime of the nebula) and the evolution at constant mass after 



the dissipation of the nebula. This phase is eventually attained 
by all planets. 

Note that the calculations of the gaseous envelope in all 
phases allows to obtain the core mass more accurately: We can 
now calculate the capture radius for planetesimals R c . dpt in the 
gas runaway accretion phase. This capture radius is necessary to 
calculate the core accretion rate Mz- In previous calculations, we 
had to extrapolate the capture radius as found in the pre-runaway 
phase into the runaway phase. The behavior of R c . dpt is discussed 
in Sect. 6.2.2. 



3.1. Structure equations 

The gas accretion rate of the planet (in the early phases before 
runaway, cf. 6.2.1) is obtained by solving the one dimensional, 
hydrostatic planetary structure equations (e.g. Bodenheimer & 
Pollack 1986). These equations are similar to those for stellar 
interiors, except that the energy release by nuclear fusion is re- 
placed by the heating by impacting planetesimals, which is the 
dominant energy source during the early formation stage. The 
other equations are the standard equations (except for the as- 
sumption of a constant luminosity within the envelope, as dis- 
cussed further down) of mass conservation, hydrostatic equilib- 
rium and energy transfer (e.g. Broeg 2009): 
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In these equations, r is the radius as measured from the plan- 
etary center, m the mass inside r (including the core mass Mz), I 
the luminosity at r, p, P, T the gas density, pressure and temper- 
ature. The quantity V(T, P) is given by 



V(T,P) 



d\nT 
dlnP 



min(V ad , V ra d) 



(3) 



where the radiative gradient (in radiative zones) in the diffusion 
approximation is given as 



kIP 



ad 



64-ncrG T 4 m ' 



(4) 



where k denotes the opacity (given by Bell & Lin 1994 and 
Freedman et al. 2008), while cr is the Stefan-Boltzmann con- 
stant. The adiabatic gradient (in convective zones) V ac j is di- 
rectly given by the equation of state SCvH (Saumon et al. 1995). 
We thus assume zero entropy gradient convection and use the 
Schwarzschild criterion to determine whether a layer is convec- 
tively unstable. 

This means that our model builds on two interlinked tradi- 
tional assumptions about the interiors of (giant) planets. These 
assumptions are mainly made for simplicity of the models, but 
are in fact not necessarily correct (Stevenson 1985; Leconte & 
Chabrier 2012): First, it is assumed, as just stated, that the in- 
teriors are adiabatic due to efficient large-scale convection. In 
this case, for Jupiter and Saturn at present time, the degree of 
super-adiabaticity (which is the fractional degree by which the 
actual temperature gradient is larger than the adiabatic gradi- 
ent) is very small (10" 8 - 10" 9 , Leconte & Chabrier 2012). 
There are however several situations where the assumption of 
an adiabatic interior may break down (Saumon & Guillot 2004). 
Additionally, at early stages of the evolution, the planets could 
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be going through regions in the surface gravity-effective tem- 
perature parameter space where the formation of H2 leads to 
a non-linear response of the atmospheric structure which influ- 
ences the degree of super-adiabaticity (Baraffe et al. 2002). We 
are currently working on including the mixing length theory (e.g. 
Kippenhahn & Weigert 1990) in order to better quantify this is- 
sue (see also Rafikov 2007). 

The second traditional assumption is that the planet consists 
of a few well separated regions, which are chemically homoge- 
neous. In our model, we assume that the envelope consists purely 
of hydrogen and helium. In reality (Stevenson 1985; Chabrier & 
Baraffe 2007), it is well possible that compositional gradients 
exist, originating for example from the dissolution of the core 
(Wilson & Militzer 2012), or the destruction of planetesimals in 
deeper layers (Mordasini et al. 2005). If a stabilizing composi- 
tional gradient is present, it is likely that the large-scale convec- 
tion breaks up into many thin chemically homogeneous convec- 
tive layers which are separated by narrow, diffusive interfaces 
with large compositional gradient (Stevenson 1979). Due to the 
large number of diffusive interfaces, the efficiency of heat trans- 
port is strongly reduced, resulting in a super-adiabatic tempera- 
ture gradient. As recently shown by Leconte & Chabrier (2012), 
such a semiconvective model is able to reproduce the observa- 
tional constraints coming from the gravitational moments and 
the atmospheric composition of Jupiter and Saturn. This shows 
that an adiabatic interior cannot be taken for granted. 

3.2. A new, simple method for the calculation of the total 
luminosity and evolutionary sequences 

In a planetary population synthesis model, the evolution of thou- 
sands of different planets is calculated, covering an extremely 
wide parameter range of planetary core and envelope masses, 
accretion rates and background nebula conditions. We therefore 
need a stable and rapid method for the numerical solution of 
these equations. We have therefore replaced the ordinary equa- 
tion for dl/dm = -TdS /dt (where S is the entropy) by the as- 
sumption that / is constant within the envelope, and that we can 
derive the total luminosity L (including solid and gas accretion, 
contraction and release of internal heat) and its temporal evolu- 
tion by total energy conservation arguments, an approach some- 
what similar to Papaloizou & Nelson (2005) and Hartmann et al. 
(1997). We first recall that -dE tot /dt - L and that in the hydro- 
static case, the total energy (neglecting rotation) is given as 



Egrav + £"int 

GM 2 



r M Gm r M 

I dm + I u 

Jo r Jm, 



dm 



2R 



(5) 
(6) 



where u is the specific internal energy, M the total mass, and R 
the total radius of the planet. The integration of the gravitational 
energy includes the core, for which we assume for simplicity a 
density which is constant as as a function of r. This is strictly 
speaking not self-consistent, as we assume that the core is dif- 
ferentiated, see Paper II. Then, the binding energy of the core 
is given as -(3/5)GM^/R core . Note that, by differentiating this 
energy with respect to time, assuming a constant density, the re- 
sulting luminosity from the growth of the core is equal to the 
accretion luminosity of planetesimals falling onto the core with 
a velocity equal to at infinity. On the other hand, the core does 
not contribute to the internal energy as we do not consider the 
thermal evolution of the core. In Eq. 6 we have introduced a pa- 
rameter which represents the distribution of mass within the 



planet and its internal energy content. This last formula is indeed 
nothing else than a definition of For example, a fully convec- 
tive, nearly isentropic star that can be approximated by a n=3/2 
polytrope would have £, = 6/7 (Hartmann et al. 1997). 

The quantity £ can be found for any given structure at time t 
with the equations above. Then one can write 



~dt^ tot ~ ^ ~ ^ M + ^ R + ^ 



£GM . (GM 2 . GM 2 . 



R 



2R 2 



2R 



(7) 
(8) 



where M = Mz + Mxy is the total accretion rate of solids and 
gas, and R is the rate of change of the total radius. This equa- 
tion corresponds to a generalization of Eq. 6 in Hartmann et al. 
(1997). All quantities except ^ can readily be calculated at time 
t. We now set 

L^C(L M +L R ). (9) 

The factor C corrects approximately for neglecting the Lg term, 
and is obtained in this way: a posteriori, one can calculate the to- 
tal energy in the new structure at t + dt, which gives the exact lu- 
minosity as L ex = -[E lot (t+dt)-E tot (t)]/dt. By setting C = L ex /L 
one obtains the correction factor so that exact energy conserva- 
tion would have occurred. As an approximation, we then use 
this C for the next time step. One finds that with this method, 
the estimated luminosity L and the actual luminosity L ex agree 
generally very well, provided that dt is small enough 1 . With this 
prescription, we can thus always calculate the total luminosity at 
t + dt, which is one of the necessary boundary conditions for the 
envelope calculations, allowing us in the end to construct evolu- 
tionary sequences. We will show below that this method leads to 
results in terms of luminosity and radius evolution which are in 
very good agreement with traditional calculations based on the 
entropy like Burrows et al. (1997) or Baraffe et al. (2003). 

In this method, the (core) luminosity due to the accretion 
of planetesimal during the formation phase, but also the release 
of energy due to the contraction of the core at constant mass - 
the dominant contribution from it during the evolution phase for 
giant planets (Baraffe et al. 2008) - is automatically included. To 
this luminosity we finally add the radiogenic luminosity L rac jio 
(see Paper II) to get the total intrinsic luminosity. Note that with 
planetary luminosity we always mean the one emitted from the 
interior, without the contribution from absorbed stellar radiation 
which makes up for Jupiter today about 40% of the total flux 
(Guillot & Gautier 2009). 

Note that we plan on further improving our scheme to esti- 
mate the luminosity at the next timestep, which is possible by 
considering the contribution from the core and the envelope sep- 
arately. 



3.3. Hydrostatic approximation 

Equation 5 assumes that the planet is always in hydrostatic equi- 
librium, in particular also during the rapid contraction phase at 
the moment the planet detaches from the nebula (see the left bot- 
tom panel in Fig. 2) which we shall call the collapse phase in- 
dependent of its true nature. To check if this assumption is self- 
consistent, we performed an order of magnitude estimate: Using 



1 The timestep dt for a calculation as in Sect. 6 is set to ~ 10 4 years 
in phase II. During the collapse phase, dt must be reduced to ~ 1 year. 
Afterwards dt can increase again, reaching ~ 1 Gyr late during the evo- 
lutionary phase. The latter two values correspond to about 0.1-1% of 
the Kelvin-Helmholtz timescale of the planet at these moments. 
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the temporal change of the total radius v(R) = dR/dt (which is of 
order -10 m/s during the fastest contraction), we have assumed 
that the contraction of the layers in the interior is homologous. 
This allows us to estimate the total kinetic energy as 

Eten = J m 2 ( v(R) r) dm - (10) 

Comparing with the total energy, we find that E^ n is always 
many orders of magnitude smaller than the total energy during 
the entire formation and evolution of a giant planet. During the 
collapse, it is smaller by about 6-8 orders of magnitude (and dur- 
ing the subsequent slow long-term contraction even by about 20 
to 30 orders). This very large difference makes one think that the 
hydrostatic approximation is justified. A caveat is that first, our 
check is only looking at global quantities and, therefore, would 
not catch if in certain layers of the planet hydrodynamic effects 
become locally important. Second, one cannot exclude a priori 
that the usage of the hydrostatic approximation itself forces the 
evolution of the planet to behave hydrostatically. 

The large security margin of 6 to 8 orders of magnitudes, 
and the fact that Bodenheimer & Pollack (1986) reached the 
same conclusion seems however indeed to indicate that the evo- 
lution does not become a dynamical collapse, but remains just 
a rapid gravitational contraction, validating our approach. We 
must however note that opposite views do exist (Wuchterl 1991). 

3.4. Simplification using dl/dr = 

The method introduced above only yields the total luminosity at 
the surface, and not l(r) which is necessary to solve the structure 
equations. We use the simplest setting, namely that dl/dr = 0. 
This seems at first sight not to be a reasonable approximation as 
it (formally) means that the complete luminosity originates in the 
core, which is, of course, not the case except for the early phases. 
One however finds that the solution of the structure equations is 
very insensitive to the specific form that / takes in the interior. 
This results can be understood due to the following physical rea- 
sons: / enters into the remaining three structure equations only 
in radiative zones of the planet (see Eq. 3 and 4). In convective 
zones, I does not enter at all. Note that this is only true for the 
assumption that convection is strictly adiabatic, which might not 
be the case, as discussed at the end of Sect. 3.1. 

However, significant radiative zones (in terms of contained 
mass) only exist during the very early phases of formation of 
the planet (see e.g. Bodenheimer & Pollack 1986). But then the 
dominant part of the luminosity is generated by the accretion of 
planetesimals by the core, so that L * L core , or in other words 
dl/dr w 0, so that our approximation is close to the exact solu- 
tion. The core luminosity due to the accretion of planetesimals 
is given as 

GM Z . 

L mK ^—^M z (11) 

^core 

where M z is the mass of the core, R core its radius, and M z is 
the accretion rate of planetesimals. This means that we use the 
sinking approximation (Pollack et al. 1996). 

During the collapse phase and during the long-term evo- 
lution when no planetesimals are accreted any more, I varies 
strongly inside the envelope, from a very small value (L ra dio) at 
the core to a typically much larger total luminosity at the sur- 
face for giant planets due to the cooling and contraction of the 
envelope. However, in this phase the planet is nearly fully con- 
vective (e.g. Bodenheimer & Pollack 1986; Guillot 2005), with 



a radiative gradient dominant over the adiabatic one by orders 
of magnitude (Alibert et al. 2005), so that the radial variation 
of / is irrelevant for the structure. A radiative layer still exists, 
but only near the very surface of the planet (e.g. Burrows et al. 
1997), containing negligible amounts of matter ( < 1%, Guillot 
& Gautier 2009). Therefore, they do not contribute significantly 
to the contractional luminosity which comes from the deeper in- 
terior (Kippenhahn & Weigert 1990), so that dl/dr = is again 
a good approximation in those parts where / matters. This also 
means that the Schwarzschild criterion, which uses I to decide if 
a layer is convective or not, remains valid. 

3.4.1 . Tests of the simplification 

Regarding the assumption that dl/dr = 0, we performed a num- 
ber of tests. We adopted different prescriptions for the form of I 
as a function of radius (or mass). First, we (arbitrarily) assumed 
a linear increase with m from l(R CO re) = ^core at the core to the 
total luminosity at the surface. 

For the second set of tests, we proceed as follows. Knowing 
the planetary internal structure at two timesteps (say t and 
t—df), we compute the Lagrangian change in the specific entropy 
S (m, t). From this, we compute Km, t) = -T(S (m, t) - S (m, t - 
dt))/dt, where T is the mean temperature for mass shell m be- 
tween t-dt and t. Then, we compute the planetary internal struc- 
ture at time t + dt by scaling this Km, t) onto the new core and 
total mass, and the new L core and total luminosity. This rescaled, 
variable l(m, t + dt) is then used to calculated the structure. 

In all cases, one finds only very small differences between 
these calculations and the simulations with a constant /. For ex- 
ample, for the time needed to go to runaway gas accretion in a 
Pollack et al. (1996) Jl like calculation (see the appendix A), 
one finds differences of 2-5% by which it takes longer in the 
dl/dr = simulation compared to the variable I case. This is a 
difference which is much smaller than those introduced by e.g. 
the uncertainties regarding the grain opacity or the planetesimal 
size. For the long-term evolution phase, the differences are vir- 
tually vanishing, because the planets are almost fully convec- 
tive. This large independence of the results on the form of I has 
also been noticed by others (A. Fortier, personal communication, 

2011) . 

Finally, we have tested the impact of the simplified lumi- 
nosity calculation during the attached and collapse phase by 
comparing with simulations based on the solution of the full 
set of structure equations using the traditional entropy method, 
obtained by a completely independent model (Broeg & Benz 

2012) . As described in details in their paper, these authors solve 
the structure equations in the quasi-hydrostatic approximation 
on a self-adaptive 1 -dimensional grid using the classical Henyey 
method (Kippenhahn and Weigert, 1990). For a simulation sim- 
ilar to the Jl case discussed in the Appendix A, we compared 
the accretion rates of gas and solids, the luminosity, and the ra- 
dius of the forming planet. Also in this case, we found very good 
agreement. 

The tests presented before may seem complicated, and one 
may wonder why we do not simply solve the whole set of equa- 
tions, replacing Eq. 2 by its standard form dl/dm = -TdS /dt . 
The reason is the following: population synthesis requires the 
computation of thousands of formation models, starting from 
different initial condition. Such a computation is done in an auto- 
matic way, using programs generating the initial conditions, then 
computing the formation model, and finally compiling the im- 
portant results (e.g. mass, semimajor axis, radius, and luminosity 
of the formed planets). In this situation, it is very important that 
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only very few simulations do not converge, since non converging 
simulations could introduce a bias in the results (e.g. if all sim- 
ulations leading to massive planets were not to converge, one 
would predict only low-mass planets). We have therefore chosen 
not to use the standard Newton-Raphson in order to compute the 
planetary internal structure, but to use a shooting method, which 
converges extremely well (less than 0.3 % of our simulation in 
the synthesis in Paper II have problems to converge using this 
method). However, the shooting method is possible in practical 
only if one knows the luminosity profile a priori, which is the 
case assuming a constant luminosity or in the two sets of tests 
presented above. 

3.5. Cold versus hot start conditions 

Finally one must take into account that a part of the released 
gravitational energy of newly accreted material is not incorpo- 
rated into the planetary structure but already lost as L acc in the 
accretion shock on the planet's surface (Bodenheimer, Hubickyj, 
& Lissauer 2000) or in the surrounding circumplanetary disk 
(Papaloizou & Nelson 2005). For matter falling in free fall onto 
the planet, the released luminosity at the surface of the planet is 
approximately given by 

GM . 

iacc = ~^~MxY (12) 
K 

where Mxy is the gas accretion rate. The same amount of energy 
is also lost in total if the matter spirals onto the planet through 
a circumplanetary disk, with the difference that some fraction 
(half of it for a Keplerian accretion disk) is already lost in the 
disk. 

We can now consider two limiting scenarios (see Spiegel & 
Burrows 2012 for intermediate "warm start" simulations). First, 
as in Papaloizou & Nelson (2005) we assume that all energy lib- 
erated in the shock is radiated away from the planet, and does not 
contribute to the luminosity inside the planet's structure. This is 
the physically likely scenario (Stahler et al. 1980; Chabrier, per- 
sonal communication, 2010). The opposite extreme is to assume 
that no energy at all is radiated away at the shock. While phys- 
ically unlikely, it is of interest to consider this case: In absence 
of radiative losses at the shock, our simulations become similar 
to the so called "hot start" calculations. This is because also for 
these "hot start" models, there is no mechanism which can act as 
a sink of entropy, which is the central role of the radiative losses 
at the accretion shock. As demonstrated Fortney et al. (2005) 
and Marley et al. (2007), the smaller initial entropy of models 
which assume radiative losses at the shock has very important 
consequences for the luminosity of young Jupiters ("cold start" 
models, see Sect. 8.2.2). 

The luminosity within the planetary structure L m for the two 
cases is given as 

_/L-L acc "Cold start" 
Unt -\L "Hot start (accreting)" UJJ 

while an observer would see the total luminosity L coming from 
both the intrinsic and accretional luminosity 2 . We have added the 

2 It is interesting to note that we should have in principle two com- 
ponents: A cooler (IR) component from the internal luminosity of the 
planet, and a harder one from the accretion shock. This is similar as 
for accreting stars. The fact that gas runaway accretion can only occur 
when there is still a sufficiently massive disk, means that the radiation 
from the planet likely gets strongly reprocessed by the surrounding disk 
material, which should make the observable signature more complex. 



suffix "accreting" for our "hot start" calculations, in order to dis- 
tinguish them from traditional "hot start" calculations (Burrows 
et al. 1997; Baraffe et al. 2003) where one starts with a fully 
formed planet (i.e. at the final mass) with some high, arbitrary 
entropy. It is clear that our treatment is a strong simplification of 
the exact physics of the accretion shock (e.g. Stahler et al. 1980). 
As already mentioned by Marley et al. (2007), three dimen- 
sional, radiation-hydrodynamic calculations resolving the shock 
physics would be necessary to get more accurate results. We will 
show in a dedicated paper (Mordasini et al. in prep.) how our 
results for the luminosities depend on a number of parameters 
which should be indicative of general trends for the behavior of 
the luminosity of newly formed giant planets. 

3.6. Possible issues for "hot" accretion anddl/dr = 

The following concern occurs when "hot" accretion is simulated 
with the dl/dr - simplification in the envelope. When the in- 
ternal energy of the infalling gas (or a part of it) is incorporated 
into the growing object rather than radiated away, this can lead 
to important modifications of the internal structure, as illustrated 
by several papers studying the effects of accretion onto low-mass 
stars (Hartmann et al. 1997; Mercer-Smith et al. 1984; Prialnik & 
Livio 1985). These studies show that under some circumstances, 
the heating of the surface layers by the rapid accretion of hot, 
high-entropy gas can reduce the temperature gradient to a degree 
that convection stops. This means that a radiative layer develops 
in the interior, even at considerable depth. 

Prialnik & Livio (1985) studied in details the reaction of an 
initially fully convective 0.2 M star to the accretion of "cold" 
and "hot" gas at a wide range of accretion rates. For the case of 
"cold" accretion, they found that independently of the accretion 
rate, the objects remain fully convective, and contract as they 
grow, indicating that the dl/dr = simplification is appropriate 
in this situation. Such a behavior of the radius is also found in 
our simulations of giant planet formation in the "cold" assump- 
tion presented in Sect. 8.2.2. If the gas keeps a fraction a/, of 
L acc (Eq. 12), different scenarios occur, depending on 07,, and 
the gas accretion rate. In this paper, we study the two limiting 
cases corresponding to 07, = 1 for completely "hot" and 07, = 
for completely "cold" accretion (see Eq. 13). 

Prialnik & Livio (1985) find that if the gas accretion rate is 
very small (M XY 5 10~ 10 M G /yr»; 3 x 10~ 5 M e /yr), the star re- 
mains fully convective also for "hot" accretion. If both the ac- 
cretion rate and 07, are high, a significant part of the entire star 
becomes radiative and the entire star expands by a large amount. 
For intermediate values of 07, and Mxy (including planetary ac- 
cretion rates expected in a protoplanetary disk < 10~ 2 M e /yr), 
~ 10 - 30% of the mass of the star becomes radiative, and the 
outer radius including the newly accreted matter expands, while 
the matter originally contained in the star contracts. An expan- 
sion of the star undergoing "hot" accretion is also seen in the 
simulations of Mercer-Smith et al. (1984). 

As discussed by Prialnik & Livio (1985, see also Hartmann 
et al. 1997), these results can be understood by comparing the 
accretion timescale of "hot" material with the thermal relax- 
ation timescale of the newly accreted material. If the accretion 
timescale is short compared to the relaxation timescale, thermal 
equilibrium cannot be established, which would be needed for 
the star to remain convective. In the "hot start" simulations (Sect. 
8.2. 1) we also see that the radii of the planets reinflate during the 
gas runaway accretion phase, after the initial collapse. These re- 
sults will be discussed in details in Mordasini et al. (in prep). It 
is clear that the presence of a deep radiative zone, which we can- 
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not catch with the dl/dr = simplification could cause strong 
departures from the results obtained with this assumption. This 
would also affect the luminosity and radius of the planets after 
formation. 

The results of Prialnik & Livio (1985) were obtained for 
spherically symmetric accretion. In reality, most of the matter 
would probably be accreted via a circumplanetary disks (e. g. 
Lubow et al. 1999), so that only a fraction of the planetary sur- 
face would be undergoing direct accretion, while the rest would 
be free to radiate. As for stars, it seems likely that this could be 
important for an accurate criterion when convection is inhibited 
by "hot" accretion (Hartmann et al. 1997). Three dimensional 
radiation hydrodynamic simulations seem necessary to clarify 
the situation. We will investigate this issue in simulations with a 
variable luminosity in future work. For the moment, this concern 
means that our results for the luminosity of young giant planets 
formed with a "hot start" (and especially at high gas accretion 
rates) must be regarded with due caution. 



found by the solution of the structure equations and is given by 
the ability of the envelope to radiate away energy so that it can 
contract (i.e. its Kelvin-Helmholtz timescale), so that new gas 
can stream in (Pollack et al. 1996). Numerically, for the shooting 
method, one iterates in this phase on the mass of the envelope. 
We use 



R 



Ra 



1 + R A /(kiiM 
t = max (p neb K neb R, 2/3) 



3rL ir 



1 int — 

l(R) = 



Sno-R 2 

Lint- 



(14) 

(15) 
(16) 



The fitting formula for the radius of the planet in Eq. 14 from 
Bodenheimer et al. (2000) reduces to Ra for Ra <k Miss^H and 
to £ii ss /?H for Miss^H «: where 



3.7. Inflation mechanism 

At the moment, we do not include any mechanism which could 
lead to the so called "radius anomaly" (Leconte et al. 2011) 
observed for many transiting Hot Jupiter. These planets have 
radii clearly larger than expected from standard internal struc- 
ture modeling as performed here. There are several possible 
physical mechanisms leading to this bloating like tidal heat- 
ing (Bodenheimer et al. 2001), dissipation of stellar irradiation 
deep in the interior (Showman & Guillot, 2002), ohmic dissi- 
pation (Batygin & Stevenson, 2010) or double diffusive convec- 
tion (Chabrier & Baraffe, 2007). Not all these mechanisms are 
yet fully understood, and we leave their study for the moment to 
dedicated works (e.g. Leconte et al. 201 1). But we note that pop- 
ulation synthesis calculations, including the effects of different 
such mechanisms, could be a fruitful way to understand which 
one reproduces best the observed radius distribution (see Paper 
II for the predicted synthetic radius distribution for a > 0.1 AU). 

3.8. Boundary conditions 

In order to solve the structure equations, boundary conditions 
must be specified, which should also provide a continuous transi- 
tion between different phases. For the formation and evolution of 
a giant planet, three different fundamental phases must be distin- 
guished: attached, detached, and evolutionary. Lower mass plan- 
ets stay in the first phase until the nebula disappears and then di- 
rectly pass into the evolutionary phase. Five boundary conditions 
are given, namely the core radius, the luminosity, the total plan- 
etary radius (in the attached phase) or the total planetary mass 
(in the detached and evolutionary phase), and the surface tem- 
perature and pressure. The four differential equations with five 
boundary condition only have a solution for a given total plane- 
tary mass (attached phase) or radius (detached and evolutionary 
phase). 

3.8.1. Attached (or nebular) phase 

At low masses, the envelope of the protoplanet is attached con- 
tinuously to the background nebula, and the conditions at the 
surface of the planet are the pressure f ne b and approximately the 
temperature r ne b in the surrounding disk. The total radius R is 
given in this regime by about the minimum of the Hill radius Ru 
and the Bondi (or accretion) radius Ra- The gas accretion rate is 



„ GM I M \ 1/3 

* A = T M^J a (17) 

The sound speed c s in the equation for the Bondi radius is cal- 
culated with the background nebula conditions. As found by hy- 
drodynamic simulations by Lissauer et al. (2009), only the inner 
part of the material inside the planet's Roche lobe is permanently 
bound to the envelope of the planet, while the upper part par- 
ticipates in the general gas flow in the protoplanetary disk. We 
therefore follow these authors in setting ka ss = 1/3. 

Regarding the temperature, Eq. 15 reflects the fact that 
the planet has to be hotter in order to radiate into the nebula 
(Papaloizou & Terquem 1999). Equation 16 finally means that 
the planet emits energy from the interior and is concurrently em- 
bedded in the protoplanetary nebula with a certain temperature, 
where r ne b is the midplane background nebula temperature as 
given by our disk model. 

3.8.2. Detached (or transition) phase 

The solution of the structure equations together with the speci- 
fied boundary conditions leads to the well known result that the 
gas accretion rate starts to increase exponentially once the core 
has grown to a mass of the order of 10 M @ (see Sect. 6). This 
means that at some moment, the gas accretion rate obtained in 
this way (i.e, through the planets Kelvin-Helmholtz timescale) 
must become higher than some external maximal possible gas 
supply. At this point, the planet enters the second phase and con- 
tracts to a radius which is much smaller than the Hill sphere 
radius. This is the "detached" regime of high-mass, runaway gas 
accretion planets. The planet adjusts its now free radius to the 
boundary conditions that are given by an accretion shock for 
matter falling at free fall velocity Vff from about the Hill sphere 
radius to the planet's surface. Probably more realistic bound- 
ary conditions would be those appropriate for the interface to 
a circumplanetary disk (Papaloizou & Nelson 2005). The maxi- 
mal possible gas accretion rate Mxy is given by how much gas 
is supplied by the disk and can flow through the gap onto the 
planet M max (e.g. Lubow, Seibert, & Artymowicz 1999). The 
calculation of this rate is specified in Section 4. The boundary 
conditions are now (cf. Bodenheimer et al. 2000; Papaloizou & 
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Nelson 2005) 



M XY = M XY ,max 



4 = 2GM 



1 1 

r~rI 



(18) 



Mxy 2p 

P = P ne b + -r-^v s + — t — max (p neh K neb R, 2/3) (19) 
4nR z 3k 



r 4 = 



3TL int 
Sno-R 2 



T 4 = (l-A)T 4 neb + Tl (20) 



The pressure in Eq. 19 has three components: the ambient pres- 
sure (which however soon after the start of the collapse becomes 
very small compared to the other terms), the ram pressure due to 
the accretion shock pv 2 , and the standard Eddington expression 
for the photospheric pressure due to the material residing above 
the t = 2/3 surface. In Eq. 20, A is the albedo assumed to have 
the same value of 0.343 as Jupiter (Guillot 2005), g is the gravi- 
tational acceleration GM/R 2 , and for the luminosity we still have 
l(R) = Li nt . In future work, we will use an albedo which depends 
on planet parameters (Cahoy et al. 2010). Numerically, one now 
iterates in this (and the next) phase on the radius R. 

For the "hot start (accreting)" models, which (artificially) as- 
sume no radiative losses at all at the shock, it is not completely 
obvious if we should not include the pressure due to the accre- 
tion shock. However, we found that the dominant term for P 
is the photospheric pressure for all accretion rates we consid- 
ered, and that the results are similar both with and without the 
ram pressure. For a gas accretion rate of 10~ 2 M ffl /yr the photo- 
spheric pressure is always a factor 15 or more larger than the 
ram pressure, while for 10 _1 M e /yr, it is still larger by a factor 
of at least about 2. This dominance of the photospheric pres- 
sure over the ram pressure stems from the low opacities at the 
typically relevant temperatures (~1500-3000 K) and densities at 
the surface of a giant planet undergoing gas runaway accretion. 
This temperature range approximately corresponds to the inter- 
val where grains have already evaporated, but where the atomic 
and molecular opacities are still low, resulting in low k ~ 0.01 
cm 2 /g even for full grain opacities. (Bell & Lin 1994; Freedman 
et al 2008). The dominance of the photospheric pressure is in- 
creased even more by the low grain opacity reduction factor 
/opa = 0.003 which is normally used in the simulations. In sim- 
ulations which are instead made with a full grain opacity, and 
M X y = 10 -1 M e /yr, the ram pressure becomes dominant by a 
factor of a few for a short time at the beginning of the detached 
phase. 

3.8.3. Evolutionary (or isolated) phase 

The last phase starts when the gaseous disk disappears so that 
the planet evolves at constant mass. During this phase, we 
use simple gray stellar photospheric boundary conditions in 
the Eddington approximation and write (Chandrasekhar 1939; 
Guillot 2005) 



P = 



3k 



1 mt 



4tto-R 2 



(21) 



T equl = 280K( I | u p (^) 7- = (1 - A< ui + 7t (22) 

and we have still 1{R) = L mt which now also equals the total 
luminosity L as L acc = (Eq. 13). This procedure means that 
our outer planetary radius R corresponds to the Rosseland mean 
t — 2/3 surface. For Eq. 22, we have assumed that the planet is 
rotating quickly and redistributing the heat from stellar irradia- 
tion over its full surface, and that the stellar luminosity scales 



as M; which applies for roughly solar-like stars on the main 
sequence. In future work, we will include boundary conditions 
which are also accurate for planets close to the parent star un- 
dergoing intense irradiation (Guillot 2010; Heng et al. 2012). 
For the moment we recall that the minimal allowed distance for 
planets in the model is 0. 1 AU where at least for giant planets, 
irradiation is not yet a very important factor. 

It is clear that the boundary conditions for the long-term evo- 
lution are described in a much simpler way than in the Burrows 
et al. (1997) or Baraffe et al. (2003) models which employ proper 
non-gray atmospheres (see Chabrier & Baraffe 1997). We have 
however found that the general agreement in terms of total lu- 
minosity or radius is very good (cf. Sections 7 and 8) and cer- 
tainly sufficient for our purpose of population synthesis. This is 
in agreement with Bodenheimer et al. (2000) who also found 
that his cooling curves using simple Eddington boundaries and 
those from Burrows et al. (1997) agree very well. 

4. Disk-limited gas accretion rate of the planet in 
the runaway regime 

The second improvement of the model presented in this first pa- 
per regards the planetary gas accretion rate in the gas runaway 
accretion phase. It occurs for supercritical cores with a mass 
larger than about 10 M e . In this phase, the gas accretion rate 
is no longer limited by the planet's envelope structure, but it is 
limited by the rate at which gas can be provided by the disk. The 
gravitational interaction of the protoplanet and the surrounding 
disk determines the accretion rate across of the growing gap (e.g. 
D'Angelo et al. 2010). 

The actual gas accretion rate Mxy must be given by the 
smaller of the two rates: 



M- 



XY 



Min(M Tk h 



M- 



XY.max I 



(23) 



where Mrkh denotes the gas accretion rate as found through the 
solution of the structure equations in the attached phase, which is 
determined by the ability of the envelope to radiate away the po- 
tential energy of the accreted matter, while MxY.max is the disk- 
limited rate. 

In our previous models, we used for the disk-limited rate 
simply 



M 



XY.max 



= M, 



disk,equi 



= ?>nvL 



(24) 



which is the viscous accretion rate in a protoplanetary disk in 
equilibrium (i.e. where the mass flux is constant as a function of 
distance from the star). This expression however neglects three 
points: (i) the presence of the planet acting as a mass sink locally 
brings the disk out of equilibrium, (ii) only a fraction of the gas 
flux through the disk will be accreted onto the planet, while some 
part will flow past it (Lubow & D'Angelo 2006). (iii) there can 
be (initially) a local reservoir of gas around the planet that can 
be accreted faster than on a viscous timescale. 



4. 1 . Non-equilibrium flux 

The inner part of a viscous accretion disk evolves in absence of a 
planet relatively quickly into equilibrium. However, the sucking 
action of the planet and photoevaporation at late stages disturbs 
the purely viscous evolution and brings the disk out of equilib- 
rium. Then, the mass flux at a distance r from the star is given 
as 

<9v£ 

M disk = 3nvT. + 67ir— (25) 
or 
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where the second non-equilibrium term becomes dominant close 
to the planet when it is in runaway gas accretion. In order to 
model the mass flux onto the planet, we have to consider differ- 
ent flow geometries of the gas in the gas feeding zone around the 
planet. We look at the outer flux Mq at Rq — R p + 1 /2R out and 
the inner flux Mi at R\ — R p - l/2R out . Here, R p is the orbital 
distance of the planet, and R out is its gas capture radius which is 
given in Eq. 14. 

4.1.1. Flow geometries 

Four different flow geometries have to be considered, where we 
count a mass flux towards the star as positive. First we can have 
a situation where both Mo and Mi are positive. This is the regu- 
lar situation. It takes place when the planet is inside the radius of 
maximum viscous couple Rmvc (or radius of velocity reversal) 
At Rmvc, the flow in the disk changes from inward accretion to 
outward spreading. Second, we can have a situation where both 
Mo and Mi are negative. This is the case if the planet is outside 
^mvc- The third case occurs when the planet is exactly at Rmvc- 
Then, decretion on both sides happens. This could in principle 
mean that only a very small amount of gas can be delivered to the 
planet. But since we take into account the local reservoir of gas 
(Sect. 4.2), this regime has no practical meaning. The fourth flow 
geometry occurs when gas flows towards the planet from both 
sides. This regime occurs due the sucking action of the proto- 
planet at the beginning of runaway gas accretion, and then again 
near the end of the disk lifetime, when the global fluxes through 
the disk are small. 

For these four geometries, the flow limited maximal gas ac- 
cretion rate MxY.max.F are calculated as 



M 



XY.max.F 



for R p < Rmvc 



Mo 

Mi for R p > Rmvc 

Abs (Mi - M ) for R p = R M vc 
Abs (Mi) + Mq for two sided accretion. 



(26) 



For the first two cases, we have assumed that the planet mi- 
grates slower than the gas flows. This is due to the fact that plan- 
ets which are in gas runaway accretion typically are also in the 
planet-dominated type II migration regime, where this condition 
is fulfilled. For the third case, we consider the net flux, and for 
the last case, the sum of both fluxes is used. 

4.2. Local reservoir 

A local reservoir of gas which is already in the vicinity of the 
planet can be accreted at a rate larger than MxY.max.Fi at least for 
some time. This is relevant at the moment when the planet just 
passes into runaway gas accretion. Such a higher accretion rate 
is important for the depth of the so-called "planetary desert" (Ida 
& Lin 2004): The higher the possible accretion rates, the lower 
the number of planets of intermediate masses (ca. 10-100 M®, 
see the discussions in Mordasini et al. 2009a and Mordasini et 
al. 2011b). 

In order to calculate the accretion rate allowed by the lo- 
cal reservoir we use an approach directly based on the local gas 
surface density, and not on the fluxes. Following D'Angelo & 
Lubow (2008) and Zhou & Lin (2007), we assume a spherical, 
homogenous, unimpeded accretion flow for a planet with a cross 
section <x cross in a gas of density p « "L/H (taken as the mean 
over the gas feeding zone) and a relative velocity v re i, 



Gas is captured at some radius R gc so that cr ctoss ~ nR gc . For sim- 
plicity, we set R gc equal to the radius defined in Eq. 14. The rel- 
ative velocity is approximately given as v re i * R gc Q. Therefore, 
the gas accretion rate in runaway, limited by the local reservoir 
MxY,max,R is approximately given as 



S , 

— QZ?L 
H gc 



(28) 



As shown by D'Angelo & Lubow (2008), such a simple descrip- 
tion of the gas accretion rate is in fairly good agreement with 
results of 3D hydrodynamic simulations. It is, however, strictly 
speaking only valid in the intermediate phase where the planet's 
Hill sphere Rh is still smaller than the disk scale height H. The 
accretion flow is no more spherical and homogenous once the 
planet has grown so massive that R H > H. Instead it has a more 
complex geometry due to the tidal interaction of the planet and 
the disk, leading to gap formation and accretion streamers (see 
Zhou & Lin 2007). However, we find that M X Y,max,R is larger 
than MxY.max.F and, therefore, the relevant quantity only during 
an intermediate phase just after the start of the disk-limited gas 
accretion phase (as illustrated in Fig. 1 below). In this phase, 
planets have typically still masses so small that Rh < H, so that 
the simple prescription is valid. 

The final accretion rate in runaway is calculated as 



M- 



XY.max 



fciubMaX (M; 



XY.max.R, ^XY.max 



x.f) 



(29) 



M- 



XY,max,R - pO~ c: 



s V re l. 



(27) 



where we have introduced a factor ki^ which takes into ac- 
count that some gas flows past the planet. Following Lubow & 
D'Angelo (2006) we use values of k Xub of 0.75 to 0.9. Finally, 
we also make sure that not more gas is accreted than the amount 
present in the feeding zone Mf ee( j, i.e. that MxY.max < A^feed/dt, 
which is however found to be automatically fulfilled with the 
above equations. 

4.3. Example for the disk-limited gas accretion rate 

To illustrate the calculation of the gas accretion rate in runaway 
described in the last section, we show in Fig. 1 an example. The 
figure shows gas accretion rates of a growing giant planet, and 
accretion rates in the disk during the phase where the planet 
starts gas runaway accretion. The planet is at about 10 AU in 
a gas and solid rich disk. The parameter ki u \, is set to 0.9, and the 
viscosity parameter a in the disk is 0.007. 

The thick solid black line shows the accretion rate of the 
planet Mxy- The thin solid black line shows the disk-limited ac- 
cretion rate, MxY.max- Shortly after 1.33 Myr, the accretion rate 
found by the solution of the structure equations Mjkh increases 
very rapidly, so that shortly afterwards, the disk-limited rate is 
reached. In this phase, MxY.max is given as £i u [,MxY,max,R- 

This 

is because MxY.max.R (shown by the magenta dashed-dotted line) 
is larger than MxY,max,F (shown by the red dotted line). We thus 
find that the local reservoir allows for some time for an accretion 
rate higher than given by the viscous rate, namely about 0.03 
M e /yr. A typical peak accretion rate of a few times 10~ 2 M ffi /yr 
has also been found in other studies (e.g. Lissauer et al. 2009). 
Later on, however, (starting at about 1.35 Myrs), we note that 
both MxY.max.R and MxY,max,F converge on approximately the 
same value. Also Mf ee( j/dt (cyan long-dashed-dotted line) takes 
then approximately the same value. 

From the curves showing the mass flux inside and out- 
side of the planet, Mi (blue dashed line) and Mo (green long- 
dashed line), we can also deduce in which flow geometry we 
are. Initially, both fluxes are positive, which means that we are 
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Fig. 1. Gas accretion rates as a function of time of a forming gi- 
ant planet, and accretion rate in the surrounding protoplanetary 
disk at the moment of transition into gas runaway accretion. The 
thick black solid line is the accretion rate of the planet Mxy- The 
thin solid black line is the disk-limited rate, Mxy.max- The red 
dotted line is MxY,max,F, the (non-equilibrium) flow limited max- 
imal accretion rate, while the rate limited by the local reservoir 
is shown by the magenta dashed-dotted line (Mxy.max.R)- The 
cyan long-dashed-dotted line shows Mf eec |/dt. The accretion rate 
in the disk inside the planet's position, M\, is the blue dashed 
line while the accretion rate in the disk outside of the planet's 
position Mo, is the green long-dashed line. 



in the standard case where the planet is inside of the radius of 
maximum viscous couple Rmvc, and the disk is globally flowing 
towards the star. At about 1.325 Myrs however, we see that M t 
becomes negative, meaning that the flow geometry changes into 
the "two sided accretion" case. This is a consequence of the ac- 
cretion onto the planet at a rate higher than the viscous one. Later 
on (at about 1 .6 Myrs, not shown in the figure) the flow geometry 
changes back into the standard geometry, and Mi is then about 
ten times smaller than Mo- This corresponds to the 10% of gas 
that are allowed to flow past the planet. Now, a new equilibrium 
has established, with the planet acting as an additional sink. 

5. Results 

We address three different subjects: First, we present a combined 
formation and evolution calculation for a Jovian mass planet at 
5.2 AU. We show the evolution of the most important observ- 
able quantities like the mass, radius and luminosity, but also the 
evolution of the central temperature and pressure. The calcula- 
tion can be compared with the calculations of, e.g., Pollack et 
al. (1996), Bodenheimer et al. (2000) or Lissauer et al. (2009). 
Second, we compare the planetary mass-radius relation for gi- 
ant planets (see also Paper II) and the influence of the core 
mass on the radius as found in our model with several other 
studies. Third, we use our model to study the luminosities of 



young Jupiters. The results for "cold starts" can be compared 
with Marley et al. (2007) or Spiegel & Burrows (2012). 



6. Example of a combined in situ formation and 
evolution simulation 

In this section, we use the extension of the planet module to 
calculate an illustrative example of a combined formation and 
evolution simulation. 

These simulations are strongly simplified in comparison to 
the calculations performed for the complete populations synthe- 
sis models (Paper II): In the full model, the moment of the onset 
of limited gas accretion (and thus the detachment from the neb- 
ula), the disk-limited runaway gas accretion rates (Sect. 4), as 
well as the final mass of the planet result in a self-consistent 
way from the evolution of the gaseous disk (Alibert et al. 2005). 
Here, the evolution of the disk is switched off. Instead, we fix 
the maximal allowed M X y to some prespecified value, and artifi- 
cially switch off accretion on a short timescale when the mass of 
the planet approaches the desired value. Also migration is com- 
pletely switched off, which is otherwise a central component of 
the formation model as it can lead to much more complex accre- 
tion histories (Mordasini et al. 2009a). But in this way we allow 
for direct comparison with numerous previous works which used 
similar assumptions as, e.g., Pollack et al. (1996), Lissauer et al. 
(2009) or Movshovitz et al. (2010). 

In Table 1, we list the most important settings used for the 
simulations. Other settings, like for example a planetesimal ra- 
dius of 100 km are, the same as in Alibert et al. (2005), except 
when mentioned. 



6.1. Grain opacity reduction factor / opa 

A new parameter in our models is the reduction factor of the 
grain opacity in the envelope relative to the ISM grain opacity, 
/opa- We have determined this factor in Mordasini et al. (2012c) 
in the following way: Movshovitz et al. (2010) coupled for the 
first time self-consistently a giant planet formation simulation 
with the concurrent calculation of the grain growth and settling 
inside the planet's gaseous envelope. They found very low re- 
sulting effective opacities, and thus very short durations of the 
(plateau) phase II (Pollack et al. 1996). Such grain evolution 
calculations are currently beyond the scope of our model, also 
because they are computationally extremely expensive. But to 
still get an approximative representation of this important effect, 
we have, for different values of / opa , and the three planetesimal 
surface densities considered in Movshovitz et al. (2010), com- 
pared the duration of the phase II as found in our model, with 
the duration found by Movshovitz et al. (2010). There is a a very 
clear relationship between / opa , the planetesimal surface density, 
and the duration of phase II. We have found that on the mean (for 
the three planetesimal surface densities), the duration of phase II 
becomes the same in both models if we reduce the interstellar 
grain opacities by a factor 0.003. 

This is a strong reduction, which is even clearly smaller than 
the previously studied "low opacity case" (Pollack et al. 1996, 
Hubickyj et al. 2005) which arbitrarily assumed a / opa = 0.02. 
Note the following two points: First, it is clear that already the 
calculations of Movshovitz et al. (2010) are simplified in several 
aspects like the usage of spherical grains in the opacity calcula- 
tions instead of aggregates, or the neglect of icy grains. Second, a 
uniform grain opacity reduction factor cannot reproduce the de- 
tailed structure found by Movshovitz et al. (2010) for the opac- 
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Table 1. Settings for the in situ formation and evolution calcula- 
tion of Jupiter. 



Quantity Value 



a[AU] 


5.2 


£ s ,o [g/cm 2 ] 


10 


MxY.max [MJyt] 


0.01 


T nc b [K] 


150 


PnA [dyn/cm 2 ] 


0.275 


Dust-to-gas ratio 


1/70 


Initial embryo mass [M e ] 


0.1 


Migration 


not included 


Disk evolution 


not included 


Planetesimal ejection 


included 


Core density 


variable 


Simulation duration 


10 10 yrs 


Grain opacity red. factor / opa 


0.003 



ity as a function of radius in the envelope. It nevertheless leads 
to envelope structures which are at least similar, especially in 
the deeper layers. These two points mean that using a uniform 
grain opacity reduction factor, calibrated via the simulations of 
Movshovitz et al. (2010), is a strong simplifications of the ac- 
tual processes in a protoplanetary envelope. It nevertheless rep- 
resents a better approximation than using arbitrarily reduced val- 
ues. The reader is referred to Mordasini et al. (2012c) for details. 

6.2. Jovian mass planet at 5.2 AU 

For this simulations, we use the same initial planetesimal surface 
density, disk pressure and temperature as Pollack et al. (1996). 
This case has been studied frequently in the literature, and we 
consider it as our nominal model. The luminosity is calculated 
as in the "cold start" assumption, i.e. we assume that the shock 
luminosity does not contribute at all to the luminosity inside the 
planetary structure. 

The formation of Jupiter is of particular importance as a 
benchmark for any planet formation model. This is because for 
no other giant planet an equally high number of detailed obser- 
vational constraints exist. In Table A.l in the appendix we list 
the most important quantities characterizing the planet like the 
total mass, core mass or luminosity at particular, important mo- 
ments in time. The data in these tables can be compared with 
similar data in previous studies mentioned above. 

The simulations presented in this section differ from those in 
Pollack et al. (1996) in a number of points, making them more 
realistic. The differences are first that we include planetesimal 
ejection, second that the core density is variable (Paper II), and 
third that / opa = 0.003, which leads to strongly reduced forma- 
tion timescale. In the Appendix A we also present simulations 
where we set these parameters to the same value as in Pollack et 
al. (1996), in order to compare quantitatively with their simula- 
tions. 

6.2.1. Formation phase: Mass 

In Fig. 2 we show the most important planetary properties as a 
function of time during the formation and the very beginning of 
the evolutionary phase once the final mass is reached. 

The top left panel shows the characteristic three stages la- 
belled with I, II and III (Pollack et al. 1996) of such classical in 
situ calculations (all three occur within the attached phase de- 



scribed in Sect. 3.8): First a rapid build-up of the core occurs 
until the isolation mass is reached, then a plateau phase follows 
characterized by a slow increase of the envelope mass which al- 
lows further core growth, and then the transition to gas runaway 
accretion is observed. 

The crossover point when the core and envelope masses 
are the same is reached relatively quickly at 0.82 Myrs due 
to the strongly reduced grain opacity. This duration is in good 
agreement with Movshovitz et al. (2010). The crossover mass 
is M CT = 16.3M e . This is in very good agreement with Pollack 
et al. (1996) with M cr =16.2M e or Hubickyj et al. (2005) who 
obtained M cl = 16.1M ffl for the same initial planetesimal surface 
density. 

Shortly after the crossover point, the gas accretion rate in- 
creases strongly because of beginning runaway accretion, so that 
it hits the maximal allowed value of 10~ 2 M e /yr at t=0.923 Myrs. 
The total mass of the planet is then 86.7 M ffi , and tkh is as short 
as ~1640 yrs. At this moment, the second detached phase begins 
and the collapse of the envelope starts, which is in fact a rapid, 
but still hydrostatic contraction as discussed in section 3.3 (see 
also Bodenheimer & Pollack 1986). This phase is indicated with 
the label D in the plot. 

At a high gas accretion rate of 0.01 M e /yr, a mass of 1M% 
is approached quickly, so that at 0.935 Myrs we start to artifi- 
cially ramp down the accretion rates (see top left panel of Fig. 
2). Shortly afterwards, at about 0.961 Myrs, the final mass of 
the planet is almost reached (99.5% of the final mass). This 
marks the beginning of the evolutionary phase, labelled with an 
E. Note that the accretion rate of planetesimals also increases 
strongly during the runaway phase because of the fast expan- 
sions of the planet's feeding zone. It reaches a maximum of 
about 7 x 10~ 4 M ffl /yr and then starts to decrease already before 
we ramp it down artificially 3 . This is due to the fact that ejection 
of planetesimals becomes increasingly important as the mass of 
the planet grows (see Alibert et al. 2005), and because the cap- 
ture radius of the planet shrinks (bottom left panel). The final 
core mass that is obtained is therefore very similar to the one that 
is found without stopping the accretion artificially, as becomes 
clear when looking at Fig. 6. 

The final total mass of the planet is 316.6 M e , while the final 
mass of the core is 32.9 M @ . In this paper, we assume that all 
accreted planetesimals can reach the core. In reality, once the 
envelope has grown above a certain mass (about 1-2 M e for 100 
km planetesimals, Mordasini et al. 2005), planetesimals cannot 
penetrate any more to the core of the planet but get destroyed 
in the envelope. Therefore, the 32.9 M e value should rather be 
associated with the total mass of heavy elements contained in 
the planet than the core mass. 

The theoretical estimates for the actual core and total heavy 
element mass in Jupiter vary. Among other factors, especially 
uncertainties in the EOS of hydrogen and helium (but also in 
the efficiency of convection) introduce considerable difficulty in 
an exact determination of the enrichment. Based on the work 
of Saumon & Guillot (2004), the total heavy element mass was 
estimated by Guillot & Gautier (2009) to lie between 10 to 42 
M e . Besides classical equations of state like the EOS SCvH 
which are based on the free energy minimization method, in the 
past few years EOS based on ab initio quantum molecular dy- 
namics calculations have become available. The results for the 



3 The point where the decreases changes from the "natural" one to the 
artificially forced one is visible in the M core line in the top right panel 
as a tiny shoulder. This is due to the tanh-like function (smoothed step 
function) we use for the artificial ramp down. 
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Fig. 2. Simulation of the in situ formation of Jupiter at 5.2 AU. Grey, vertical lines show the major phases, labelled in the top left 
panel: I, II, III during the attached phase, D for the detached phase, and E for the evolutionary phase. The top left panel shows the 
evolution of the core mass (red solid line), the envelope mass (green dotted line) and the total mass (blue dashed line). The top right 
panel shows the accretion rate of solids Mz (red solid line) and of gas Mxy (green dotted line). The limiting gas accretion rate is 
fixed to 10~ 2 M e /yr. The bottom left panel shows the evolution of the core radius R cole (red solid line), the total radius R (blue dashed 
line) and the capture radius for planetesimals R c . dpt (green dotted line). The bottom right panel shows the luminosity of the planet in 
present day intrinsic luminosities of Jupiter (Z/i,_=8.7xlO~ lo L ). The red solid line is the total luminosity L, the blue dashed line is 
the internal luminosity L mt and the green dotted line is the core luminosity L core . 



core and total heavy element mass of Jupiter derived from this 
method, however, vary: Militzer et al. (2008) find a massive core 
of M core = 16±2M ffi and M z = 20+4M e , which is quite different 
from the result of Nettelmann et al. (2012) of M core = - 8M e 
and M z = 28 - 32M ffi . Finally, Leconte & Chabrier (2012) who 
use again EOS SCvH, but allow for semiconvection due to stabi- 
lizing compositional gradients, find high M z = 41 -63M ffl as the 
interiors are hotter. These results demonstrate that the investiga- 
tions regarding the composition of Jupiter have not yet settled on 
a definitive conclusion. 



In our model, the accretion of the core occurs in two steps: 
about 10-15 M @ are accreted in phase I and phase II, while the 
rest is accreted in the gas runaway accretion phase. Then, the 
mass of the envelope is already high, and planetesimals should 
be destroyed in the envelope. Therefore, we can estimate that the 
actual core mass should rather be of order 12 M m (this is when 
the envelope mass reaches 2 M m ), while the remaining 21 M @ of 
metals should be dissolved in the envelope, unless they sink to 
the core-envelope interface (Pollack et al. 1996) . 
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Fig. 3. Left panel: Radius R (blue dashed line), core radius R cole (red solid line) and capture radius R capl (green dotted line) as a 
function of time. The inset figure is a zoom-in onto the late evolution and shows the radius as found in this work (solid line), in 
Baraffe et al. (2003) (dashed) and Burrows et al. (1997) (dotted). Right panel: Total luminosity L as a function of time (red solid 
line). The dashed line shows Baraffe et al. (2003) and the dotted one is Burrows et al. (1997). 



6.2.2. Formation phase: Radius 

The bottom left panel of Fig. 2 shows the evolution of the total 
radius, the core radius and the capture radius. The total radius 
is initially (during the attached stage) very large, as it is approx- 
imately equal to one third of the Hill sphere radius. When the 
limiting gas accretion rate is hit, the planet detaches from the 
nebula and collapses to a radius of initially about 2 to 3 Jovian 
radii. The shape of the curve showing the total radius is very 
similar to that found by Lissauer et al. (2009). 

The capture radius (which is found by integrating the orbits 
of planetesimals inside the envelope, see Alibert et al. 2005) is 
during phase II much larger than the core radius, namely by a 
factor between ~10 and 60. This is very important for the core 
growth, as the capture radius enters quadratically into the growth 
rate. After the envelope has started to collapse, /? capt becomes 
equal or slightly less than the total radius. This is due to the fact 
that in the collapse the envelope structure becomes much more 
compact. Before the collapse, it is characterized by very large 
scale heights of thousands to several ten thousands of kilometers, 
where planetesimals get captured in a soft catch. Durning the 
collapse the planet develops a well defined surface (i.e. a much 
smaller scale height near the surface) at which planetesimals are 
captured. 

Looking at the core radius, we note that in the gas runaway 
accretion phase when the planet grows very quickly in mass, it 
first increases a bit in size, but then decreases, even though the 
core grows in this phase by about 20 M m . The reason is that 
the external pressure on the core by the accreted gas increases 
so much that the resulting strong increase of the core density 
over-compensates the mass increase (see Paper II for a dedicated 
discussion). 

The bottom right panel of Fig. 2 shows the core lumi- 
nosity L core , the total luminosity L and the luminosity within 
the planetary structure L mt . The difference L - L lnt is thus the 



shock/accretion luminosity L acc . The first peak in the curve is 
due to the rapid accretion of the core, and the second to the com- 
bined effects of runaway gas accretion and envelope collapse. 
During the second peak, for about 2 x 10 4 years, the accretional 
luminosity is larger than L; nt by up to a factor ~3. For more mas- 
sive planets discussed below, this ratio can become much larger. 
The second luminosity maximum occurs at 0.94 Myrs with 
L = 2.6xl0 6 L% corresponding to log(L/L ) = -2.65. Lissauer 
et al. (2009) find in their simulations with the same limiting gas 
accretion rate peak luminosities of about log(L/L Q ) = -2.35, 
which is similar to our value. 

We note that the classical Pollack et al. (1996) picture with 
three well defined phases might not be a realistic formation sce- 
nario for Jupiter. Several factors including migration (Alibert et 
al. 2004), a slower core growth (Fortier et al. 2007) or the open- 
ing of a gap in the planetesimal disk (Shiraishi & Ida 2008) 
might suppress in particular phase II. We had two reasons for 
investigating the classical scenario: First, we want to compare 
our model with previous studies and second, we are primarily 
interested in the collapse and evolutionary phase which should 
be similar in modified scenarios. 

6.2.3. Evolution phase: Radius and luminosity 

Figure 3 shows the radius and luminosity of the planet as a 
function of time, now including the long-term evolution over 
Gigayears when the mass is constant. The evolution now oc- 
curs slowly by gradual contraction and cooling. For compari- 
son, the temporal evolution of the radius and the luminosity of 
a 1 M% planet as found by Baraffe et al. (2003) and Burrows et 
al. (1997) is also shown. The temporal zero point in these mod- 
els has been associated with the time the planet reaches its final 
mass in our calculations. Baraffe et al. (2003) and Burrows et 
al. (1997) use non-gray atmospheric boundary conditions, and 
the standard method based on the -TdS jdt term to calculate the 
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Fig. 4. Left panel: Central conditions as a function of time. The red solid line is the temperature at the core-envelope interface, 
while the dotted black line shows the pressure at this point. These two curves belong to the left axis. The green long-dashed and 
the blue dashed line are the mean core and the central gas density, respectively. They belong to the right axis. Right panel: Surface 
conditions as a function of time. The red solid line is the temperature at the surface, T sm f, while the dashed-dotted line visible at late 
times shows the temperature due to the intrinsic luminosity only. The black dotted line is the background temperature. These lines 
belong to the left axis. The blue dashed line is the (total) pressure P. The green long-dashed line shows the contribution from the 
ram pressure during the gas runaway accretion phase. These lines belong to the right axis. 



evolution. As we can seen, there is good agreement among the 
models. The differences between our model and the other two 
ones are of the same order as the mutual differences between 
Baraffe et al. (2003) and Burrows et al. (1997). Note that the 
different internal compositions (mass of heavy elements) have 
a certain influence on the exact value of R, so that we cannot 
expect exactly identical results. 

The planet has the following properties after 4.6 Gyrs: A 
radius of 0.99R% and a luminosity of 1.1 3L% corresponding to 
log(L/L ) = -9.01. At the same age, Burrows et al. (1997) find a 
radius of 1.11R% and a luminosity of \.52L%. Larger values are 
expected, as the models in Burrows et al. (1997) are core-less 
models. 



6.2.4. Central conditions 

In Figure 4, left panel we show various central conditions dur- 
ing both the formation and evolution phase. With "central" con- 
ditions we mean here the envelope-core boundary, and not the 
actual center of the planet. 

The solid red line in the left panel shows the central temper- 
ature. It belongs to the scale on the left. Initially, the temperature 
at the bottom of the envelope is about 3000 K, rising to values 
between 10000 and 20000 K during phase II. When gas run- 
away accretion starts, and the envelope collapses, there is a sharp 
upturn of r cen t. The maximal value of 69 700 K is reached at 0.98 
Myrs, afterward it declines again. At an age of 4.6 Gyr, r cent has 
fallen to about 17 600 K, which is in good agreement with the 
estimates of Guillot (1999) ranging from 15 000 to 21 000 K. 

The second quantity for which the scale is shown on the left 
is the central pressure, indicated by the black dotted line. It is 
found to have a value of 43 17 GPa at the current age of the Solar 
System. Guillot (1999) gives a value of 4000 GPa. 



The other two lines in the figure belong to the y-axis on the 
right. The green long-dashed line is the mean density of the core, 
while the blue dashed line is the gas density at the envelope-core 
boundary. As described in Paper II, the core consists of iron, 
followed by a silicate layer and finally an ice layer, with an ice 
fraction expected for a body forming outside the iceline. During 
phase II, p core lies between 3 and 4 g/cm 3 . This is close to the 
typically assumed constant value of 3.2 g/cm 3 in other formation 
calculations (e.g. Movshovitz et al. 2010). During the collapse, 
the external pressure exerted onto the core increases dramati- 
cally, so that the the core density increases to about 10 g/cm 3 , 
leading to the mentioned shrinking of the core size, despite its 
growth in mass. At 4.6 Gyr, we find p core = 14.3 g/cm 3 . The 
central gas density quantitatively follows a similar pattern, start- 
ing of course at a smaller value of about 0. 1 g/cm 3 during phase 
II. At 4.6 Gyr is has reached a fluid like density of about 4 g/cm 3 . 

6.2.5. Surface conditions 

The right panel of Fig. 4 shows the conditions at the surface 
of the planet. The solid red line is the surface temperature T sar f 
calculated as 

r - rf = r4 + 4^' (30) 

where T is defined in Eq. 20. It thus contains the contributions 
from the intrinsic luminosity, the accretional luminosity which 
is however important only for a short period (see Fig. 2 bottom 
right panel) and the absorbed stellar light. When the accretional 
luminosity is negligible, r sur f = T, which is the temperature used 
as the outer boundary condition for the structure. When both L[ nt 
and L acc are negligible, T sur f = (1 - A)r equ j. This happens for 
a Jovian planet at 5.2 AU after > 10 Gyrs when the intrinsic 
luminosity of the planet becomes small compared to the inci- 
dent stellar flux. The background temperature is indicated by 
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the black dotted line. The background temperature during the 
attached and detached phase is equal to the (constant) nebula 
temperature r ne b = 150 K, while in the evolutionary phase, it 
is given by the equilibrium temperature with stellar irradiation 
Jequi at 5.2 AU (Eq. 22). The two temperatures belong to the 
left scale. One finds that in phase I, r slll f can become about a 
factor 2.5 larger than T as b, as L is quite large in this phase. In 
phase II, T sm i as T ne b, in agreement with Papaloizou & Nelson 
(2005). During the collapse, there is a strong increase of r sur f, 
which reaches a maximal value of 2050 K. At the moment the 
final mass is reached, T sur f « 1350 K. At the current age of 
Jupiter, T sm f = 126 K, while the measured value is about 124 
K (Guillot & Gautier 2009). We see that at late times, the de- 
crease of r sur f starts to flatten out as it approaches the constant 
value of (1 - A)r equ i. In the figure the dashed-dotted red line 
is (L ml l{Ana-R 2 )) l l A i.e. the temperature due to the intrinsic flux 
only. This quantity continues to fall. Note that the nomenclature 
for the different temperatures is unfortunately not homogenous 
in the literature: Our quantity r sur f corresponds during the evo- 
lutionary phase to 7\herm according to the definition in Baraffe et 
al. (2003), but to r eff in Fortney & Nettelmann (2010). 

In addition, we show the pressure P at the surface of the 
planet. This curve belongs to the right axis. During the attached 
phase, it is simply equal to the (constant) nebula pressure (Eq. 
14). Then, at the moment that the planet detaches, P increases as 
^edd = 2/3 g/K becomes large. The dynamic/ram pressure pv 2 
also contributes to P in this phase (Eq. 19). The ram pressure 
remains however always 1-2 orders of magnitude smaller than 
Pedd which is related to the low grain opacity (Sect. 3.8.2). At 
late times, after the accretion is over, we see that P has a pecu- 
liar shape which is given by opacity transitions. At 4.6 Gyr, the 
pressure is about 0.18 bars. 

6.2.6. Atmospheric pressure-temperature profile 

In Figure 5 we show the pressure as a function of temperature 
for the uppermost part of the planet down to a pressure of 1 kbar 
at different moments in time. The evolutionary sequence starts at 
the first structure (at the right) at a time of 1 Myr (i.e. shortly af- 
ter the planet has reached the final mass), while the last structure 
(at the left) corresponds to 4.6 Gyr. Structures are shown at ir- 
regular intervals between these times. We recall that the surface 
of the planet and thus the upper end of the lines corresponds to 
the r = 2/3 surface. During the evolutionary sequence, the ra- 
dius of the planet shrinks from 2.36 to 0.99 R%, which causes an 
increase of the gravitational acceleration g from 444 at the be- 
ginning to 2546 cm/s 2 at 4.6 Gyrs. In the figure, the thick black 
lines indicate radiative zones, while the blue dotted line shows 
the p — T structure measured by the Galileo probe (Lodders & 
Fegley 1998). We find a very good agreement. 

The figure can be approximately compared with Fig. 2 in 
Burrows et al. (1997). These authors use non-gray atmospheres, 
but indicate the location which corresponds to our outer sur- 
face. One must also take into account that they keep the gravity 
at a constant value of 2200 cm/s 2 so that we cannot expect to 
see identical results, especially at early times. But a comparison 
shows that the two p — T profiles nevertheless have a similar gen- 
eral structure. In both works, there exists for example a detached 
radiative zone at temperatures around 1500 K. 

In summary we see from the in situ formation and evolution 
calculation for Jupiter that our combined formation and evolu- 
tion model, despite the simplifications, leads to a planet with 
basic properties in very good agreement with the most impor- 
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Fig. 5. Evolutionary sequence of the pressure-temperature pro- 
files near the surface of the planet. The first profile (on the right) 
corresponds to t — 1 Myr, while the last one (on the left) is at 
t = 4.6 Gyrs. The thick black lines show radiative zones. The 
blue dotted line is the profile measured by the Galileo space 
probe (from Lodders & Fegley 1998). 



tant observational constraints coming from Jupiter, and that the 
new simple method to calculate evolutionary sequences leads to 
equivalent results as the traditional entropy based method. 

7. Radii 

We now generalize the results from the last chapter concerning 
the radius and study the radii of giant planets of different masses 
and of different compositions. The goal is to validate our model 
by comparison with existing work (e.g. Fortney et al. 2007 or 
Baraffe et al. 2008). 

We do this by performing the same combined formation and 
evolution calculations as shown for Jupiter, with the only differ- 
ence that we vary the initial planetesimal surface density (which 
will eventually lead to different core masses) and the moment 
when we terminate the gas accretion, so that we get different 
total masses. Otherwise, the calculations are identical, which 
means for example that our calculations apply for planets at a 
distance of 5.2 AU from a solar like star. 

Such calculations are shown in Fig. 6. The plot shows for 
planets of final masses of 0.15, 1, 2, 5 and 10 M% the total and 
core mass as a function of time. The lines for the 1 Jupiter mass 
planet are the same as in Sect. 6. The gas accretion rate in run- 
away is 0.01 M ffl /yr in all cases. Therefore, more massive planets 
reach their final mass later. It is interesting to note that all plan- 
ets have a final core mass of about 33 M®, except for the lowest 
mass planet (total mass 0.14M% = 44.5M ffi ) which has a core of 
about 18 Earth masses. The lower core mass in this case is due 
to the fact that for this planet, gas and solid accretion are (ex- 
ternally) ramped down before it reaches the gas runaway phase. 
The nearly identical core mass of all other planets is in contrast 
not externally imposed. It is rather a natural consequence of the 
decrease of the capture radius at the moment when the planet 
collapses (which happens always at the same mass, independent 
of the final mass) and the ejection of planetesimals which be- 
comes important as the planet grows in mass (Sect. 6.2.1). 
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Fig. 6. Total mass (solid lines) and core mass (dotted lines) as 
a function of time for planets with final masses of 0.14, 1, 2, 5 
and 10 M%. For the lowest mass planet we shut down accretion 
before it passes into the runaway gas accretion regime. Note that 
the four more massive planets have a nearly identical core mass, 
which is not externally imposed, but a natural consequence of the 
decrease of the capture radius and the increase of planetesimal 
ejection. 



4.5 

^ 4 
H 



m 
o 



3.5 



2.5 



1 1 


1 1 1 


■ 

/20Mj 






//lOMj 


- 




5 Mj 




Ay/ 1 Mj 




- 


y^yyy o.i4Mj 


- 




■ i i 





2 4 6 

log(P/bar) 



10 



12 



Fig. 7. Internal temperature-pressure profiles for planets with 
different masses as labeled in the figure, at 5.2 AU from the Sun, 
and an age between 4 and 5 Gyrs. The core masses of the plan- 
ets are about 33 Earth masses, except for the Q.14M% (44.5 M e ) 
planet which has a core of about 18 M m . 



Figure 7 shows the internal pressure-temperature structure 
of these planets (plus also of a 20 M% planet) at an age of 4 to 

5 Gyrs. The left end of the lines corresponds to the surface of 
the planet, while the right end corresponds to the envelope-core 
interface. The 0.14 M% planet has near the surface a significant 
radiative zone. Otherwise the planets are nearly fully convective, 
and characterized by a single adiabat. Near a pressure of about 1 
Mbar we see a change in slope which comes from the molecular 
to metallic transition of hydrogen, also visible in a similar figure 
in Guillot & Gautier (2009). 

7.1. Mass-radius relation 

The mass-radius relation of (giant) planets has been studied for 
a long time (e.g. Zapolsky & Salpeter 1969). The interest in the 
relation lies in its connection to the composition of the planet and 
the state of matter in its interior. For recent reviews, see Chabrier 
et al. (2009) or Fortney et al. (2010). 

The general result qualitatively already found by Zapolsky 

6 Salpeter (1969) is that the M-R relationship in the giant planet 
regime is characterized by a local maximum in R. This behav- 
ior can be understand with polytropic models with a polytropic 
index that increases with mass. This change is in turn due to 
the increasing importance of degeneracy pressure of electrons 
relative to the classical coulomb contribution of ions with planet 
mass (e.g. Chabrier et al. 2009). Another general result is that the 
radius of giant planets decreases with core mass and increases 
with proximity to the star (e.g. Fortney et al. 2007; Baraffe et al. 
2008). 

In Figure 8 we show the mass-radius relation for planets with 
masses between 0.14 and 20 M\. This are the same models as 
in Fig. 6 and 7. Except for the lowest mass planet with a core 



of about 18 M @ , the core mass is approximately 33 M @ . The 
radii are shown at 0.1, 1, 4.6 and 10 Gyrs. The planets were 
calculated using the cold start assumption, which however plays 
only a role at t = 0.1 Gyrs and for planets more massive than 
5 M%. In the top right corner of the figure we show with the 
magenta dotted line the radius for simulations performed with 
the "hot start (accreting)" scenario. The radii are as expected 
larger, by about 0.07 R% for the most massive planet. For all other 
cases the imprint of the formation is lost by ~0. 1 Gyrs, and the 
radii are virtually identical. 

The radii decrease with time, as expected. At all times, there 
is a maximum of the radius at a mass of about 3-5 M%, in agree- 
ment with e.g. Chabrier et al. (2009) or Fortney et al. (2010). 
For comparison, we also show the results from a more detailed 
model by Fortney et al. (2007). The very light blue dotted lines 
show the radius as found by them at an age of 4.5 Gyrs for plan- 
ets with a core of 25 M e (and thus similar as here) at orbital 
distances of 1 AU (upper line) and 9.5 AU (lower line). We see 
that these lines bracket our red line at almost all masses, which 
is exactly what we expect when the two models yield very sim- 
ilar results. Only for the most massive planets > 10M^ we see 
that our models yields slightly larger radii, by a few 0.01 R%. 
Also the lowest mass planet has in our simulations a somewhat 
larger radius than the counterpart in Fortney et al. (2007). This 
is however not a surprise as it has a smaller core by about 6 M @ , 
which is relevant due to the small total mass of about 44 M m (cf. 
Fortney et al. 2007). 

7.2. Influence of the core mass 

The fact that a more massive core for a given total mass leads 
to a smaller radius has important implications for the study of 
planets. It allows us for example (within certain limits) to infer 
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Fig. 8. Radius of giant planets as a function of mass, for differ- 
ent moments in time. The lines correspond from top to bottom to 
t-Q.l (magenta, long-dashed), 1 (blue, dashed-dotted), 4.6 (red, 
solid) and 10 Gyrs (green, short-dashed). All planets are at 5.2 
AU from the sun, and have a core of 33M ffi except for the 0.14 
M% planet which has a core of about 18 M e . The two dotted, 
very light blue lines show the result from Fortney et al. (2007) 
for planets with cores equal 25 M @ at 1 AU (upper line) and at 
9.5 AU (lower line), at 4.5 Gyrs. The dotted branch of the ma- 
genta line in the top right corner shows the result for "hot start 
(accreting)" simulations. Otherwise the radii are virtually iden- 
tical for the "hot" and "cold start" scenario at these late times. 



the bulk composition of transiting extrasolar planets, which in 
turn constrains formation models (see Paper II). It is for example 
clear that giant planets formed by the core accretion mechanism 
must be enriched in heavy elements like it is also the case for 
Jupiter. For the competing formation model, the disk instability 
model, the situation is in contrast less clear (Boley et al. 2011). 

Another important application is to connect the inferred core 
masses of transiting extrasolar planet with their formation en- 
vironment. In this context, it was found that the core mass of 
giant planets and the metallicity of their host star is positively 
correlated (Guillot et al. 2006; Burrows et al. 2007). This is re- 
produced by planet population synthesis calculations based on 
the core accretion paradigm (Mordasini et al. 2009b). Recent 
findings (Miller & Fortney 201 1) even indicate that all transiting 
giant planets contain a minimum amount of 10-15 M ffl of metals, 
as predicted earlier by the core accretion mechanism (Mordasini 
etal. 2011c). 

It is therefore important to study the dependence of the to- 
tal radius on the core mass and to compare this with previous 
studies. Figure 9 shows the radius as a function of core mass for 
planets with a total mass of 1 M\. The red solid line shows the 
result from this work for five planets at an age of 4.6 Gyrs. We 
see that for an increase of the core mass from ~ 10 to ~ 50 M e , 
there is a reduction of the radius by about 0.07 R\, We also show 
the results from Fortney et al. (2007) and Baraffe et al. (2008). 
Our result for a — 5.2 AU are located between the curves for 
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Fig. 9. Radius of a 1 Jupiter mass planet as a function of the core 
mass for different models. Red solid line: this work, for a = 5.2 
AU and t = 4.6 Gyrs. Blue dashed line: Fortney et al. (2007), 
a — 1 AU and t = 4.5 Gyrs. Green dotted line: Fortney et al. 
(2007), a = 9.5 AU and t = 4.5 Gyrs. The magenta long-dashed 
line is from Baraffe et al. (2008), for an isolated object at an age 
of 5 Gyrs. 



1 and 9.5 AU from Fortney et al. (2007) at 4.5 Gyrs, with very 
similar slopes of the curves. From this we conclude that the two 
models yield very similar results. 

An interesting aspect of the combined formation and evolu- 
tion calculations is that we can directly relate the radius (respec- 
tively the core mass) with the initial surface density of planetesi- 
mal during formation, which is not possible for purely evolution- 
ary models. In the current case, the five different core masses cor- 
respond to planetesimal surface densities of 3.5, 5, 7.5, 10 and 
15 g/cm 2 . We must however take into account that there is no 
unique mapping between planetesimal surface density and core 
mass, as other factors like migration or the distance where the 
planet forms also influence this quantity. 

7.3. Comparison and limitations 

We conclude from the different comparisons that our upgraded 
model allows to get evolutionary sequences and thus radii with 
very good accuracy starting with self-consistent initial condi- 
tions. 

It is however also clear that we have mostly addressed rela- 
tively "simple" cases of massive, gas-dominated planets at large 
distances from the star. Already in the solar systems, there are 
planets which could not be very accurately modeled with a sim- 
ple model as used here: It is for example well known that for 
Saturn, He demixing and settling is necessary to explain the 
measured luminosity (e.g. Hubbard et al. 1999). We must also 
expect that for low-mass, core-dominated planets, possibly with 
strongly enriched envelopes, we would get less accurate cool- 
ing curves and radii, as we do not model the thermodynamics of 
the core, assume that all solids are in the core (cf. Baraffe et al. 
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2008), and also ignore effects of the composition on the opac- 
ity. Other limitations apply for strongly irradiated planets. These 
effects will be added in later modifications of our model. 



8. Luminosities 

The other fundamental observable quantity apart from the ra- 
dius which we obtain from our upgraded model is the plane- 
tary luminosity. It is obviously central for direct imaging obser- 
vations. For transits observations, the discovery of giant plan- 
ets with radii much larger than expected from standard internal 
structure modeling was a surprise. Regarding the direct imag- 
ing technique, there is on theoretical grounds an open question, 
too: what is the luminosity of young Jupiters? Currently, two 
classes of models are discussed in the literature: "hot start" mod- 
els (Burrows et al. 1997; Baraffe et al. 2003), and "cold start" 
models" (Fortney et al. 2005, Marley et al. 2007). In the "hot 
start" model, one assumes as the initial state an already fully 
formed planet at an (arbitrarily) high entropy state, and thus a 
large radius and luminosity. In the "cold start" scenario one as- 
sumes that the planet is gradually built up by accretion of gas 
through an accretion shock on the surface of the planet. The ra- 
diative losses of the liberated gravitational energy at the shock 
lead to a lower entropy, and thus lower luminosity and radius. 

8.1. Relationship of cold/hot start and core 
accretion/gravitational instability 

The importance of the two different initial states stems from the 
fact that the luminosity of young planets is an observable quan- 
tity establishing links to the physical mechanisms occurring dur- 
ing the formation. The currently known directly imaged planets 
seem to be incompatible (Janson et al. 201 1; Marley et al. 2012) 
with the "cold start" model of Marley et al. (2007), since the ob- 
served temperatures and luminosities are higher than predicted 
by this model. It should be kept in mind that the model of Marley 
et al. (2007) studies the limiting case of a shock that is 100% ef- 
ficient in radiating away the accretional energy (and makes a 
number of assumptions concerning the formation process, cf. 
below). This limiting case is studied as detailed, possibly multi- 
dimensional radiation-hydrodynamic calculations of the shock 
properties have not yet been conducted. Recently, Spiegel & 
Burrows (2012) have therefore studied a broad range of inter- 
mediate "warm-start" models showing that this leads to signifi- 
cant differences in the brightness of the planets depending on the 
mass and spectral band. 

The related question about the initial thermodynamic state 
of newly formed brown dwarfs and stellar objects has been stud- 
ied in various works (e.g. Prialnik & Livio 1985, Hartmann et 
al. 1997; Baraffe et al. 2009; Commercon et al. 2011). These 
studies indicate that depending on the fraction of accretion lumi- 
nosity radiated away or absorbed by the protostar, the formation 
through gravitational collapse can lead to bright, extended ob- 
jects ("hot accretion" leading to a "hot start") or faint and com- 
pact objects ("cold accretion" leading to a "cold start"). The lat- 
ter is the case if the accretion shock during the formation process 
is supercritical, i.e. if all the accretion shock energy is radiated 
away, as found recently by Commercon et al. (2011) for accre- 
tion on the first Larson core. 

For giant planets, two formation mechanisms are currently 
discussed. If a "cold" or "hot start" can preferentially be asso- 
ciated with one of the two formation mechanisms, then the un- 
derstanding of the luminosity of young Jupiters has important 



implications for the understanding of giant planet formation as a 
whole. 

8.1.1. Core accretion 

The first proposed formation mechanism is core accretion, which 
is the paradigm underlying this work. A number of observational 
findings indicate that core accretion represents the dominant for- 
mation mechanism at least for planets inside a few AU (e.g. 
Mordasini et al. 2012a). Whether core accretion can also explain 
the directly imaged, massive planets at orbital distances » 10 
AU is the subject of an ongoing debate (e.g. Dodson-Robinson 
et al. 2009; Kratter et al. 2010). The core accretion, "cold ac- 
cretion" simulations of Marley et al. (2007) lead as mentioned 
to very low initial luminosities, incompatible with the observed 
values (Janson et al. 2011; Marley et al. 2012). We will how- 
ever show in Mordasini et al. (in prep.) that this result (found for 
a specific set of parameters during the formation process) does 
not mean that planets formed by core accretion in a "cold accre- 
tion" scenario must necessarily (for all reasonable parameters) 
have luminosities as low as in Marley et al. (2007) (see also the 
recent calculations of Spiegel & Burrows 2012). Furthermore, 
the results in the next two sections show that the formation of 
giant planets via core accretion, but with "hot accretion" leads 
to high initial luminosities which very quickly converge on the 
luminosity of classical "hot start" calculations (Burrows et al. 
1997; Baraffe et al. 2003). We thus see that depending on cur- 
rently ill constrained physical mechanisms occurring during the 
formation phase (in particular the nature of the accretion shock), 
the luminosity of planets formed by core accretion can probably 
vary significantly. 

8.1.2. Gravitational instability 

For the second formation mechanism, the gravitational insta- 
bility (GI) model, the question about the associated luminosity 
might also be relatively complex as indicated by the explorative 
discussion in this section. More solid results are expected once 
multidimensional, radiation-hydrodynamic calculations follow- 
ing the entire formation by gravitational instability from the first 
instability in the disk to a compact planet at the final mass will 
become feasible. 

In a qualitative way, we could expected that the post- 
formation luminosity of planets formed by GI depends on the 
specific way these objects grow in mass: Hydrodynamic simu- 
lations show that the initial mass of the clump which becomes 
self-gravitationally unstable is typically 3-20 M% (e.g. Forgan 
& Rice 2011; Helled & Bodenheimer 2011; Zhu et al. 2012). 
These masses roughly correspond to the ones expected analyt- 
ically from the mass contained in the disk patch with a length 
scale similar to the most unstable mode. The objects that form at 
this stage are self-gravitating, relatively cold, extended objects. 
Their radius is comparable to the Hill sphere radius (e.g. Zhu 
et al. 2012). For example, a 10 M% clump at 50 AU from a 1 
M star which may represent a typical outcome (e.g. Janson et 
al. 201 1) has a Hill sphere radius of about 7.4 AU. These clumps 
consist of molecular hydrogen and, once formed, contract slowly 
(e.g. Helled & Bodenheimer 2011). In their nature they corre- 
spond to the first Larson core of star formation. 

The mass contained in such an object has not gone through 
an accretion shock, since they form from the spontaneous self- 
gravitational binding of a extended part of the disk, typically 
within an over-dense spiral arm (e.g. Boley et al. 2010). We 
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therefore expect that the specific entropy of the gas in the clump 
should be high. This is confirmed by numerical simulations: T. 
Hayfield (personal communication, 2012, see also Galvagni et 
al. 2012) finds that the specific entropy in a clump taken from 
the radiation-hydrodynamic simulations of Boley et al. (2010) is 
about 15 ^B/baryon at the moment when the clump should un- 
dergo the second collapse (cf. below). If this is representative 
of the entropy of the final object, this corresponds to a "hottest 
start" (Spiegel & Burrows 2012). This allows to make a first sup- 
position: If the final mass of a GI planet is the same (or similar) 
as the mass of the initially unstable clump (an approximation 
which is regularly made, e.g. Janson et al. 2011), then we may 
expect that these objects have a "hot start", simply because there 
is no accretion shock involved during their formation. 

However, since the clump necessarily forms in a massive 
disk, it appears more likely that the clump continues to accrete 
mass (in principle, gap formation can stop accretion, e.g. Zhu et 
al. 2012). In this case, we have to distinguish two sub-scenarios. 
As the first core contracts quasi-statically, its interior heats up, 
until a central temperature of about 2000 K is reached. Then, 
molecular hydrogen dissociates, and a dynamical collapse of the 
entire clump occurs (e.g. Bodenheimer et al. 1980). This cor- 
respond to the formation of the second core in star formation. 
The timescale until this second collapse occurs is estimated to 
be roughly 10 3 to 10 4 years, with more massive clumps collaps- 
ing earlier (Helled & Bodenheimer 201 1). 

In the first sub-scenario the clump grows to its final mass 
predominately before it collapses. Accretion in this phase is effi- 
cient, since the cross section is approximately given by the large 
Hill sphere radius. This results in high to very high accretion 
rates of 10" 7 to 10- 4 M o /yr (Boley et al. 2010; Zhu et al. 2012). 
The object therefore may grow quickly out of the planetary mass 
regime (Kratter et al. 2010). Independently of this issue, for the 
"hot" vs. "cold accretion" question, it is relevant in which way 
the gas is incorporated into the clump, in particular if it is ac- 
creted via a supercritical shock. The escape speed from the Hill 
sphere of a 10 M% object at 50 AU is about 1.5 km/s. For com- 
parison, in the calculation of the accretion onto the first (stellar) 
core in Commercon et al. (201 1), the shock velocity is about 2.4 
km/s, which is comparable. Because we are in a shearing disk, 
both Boley et al. (2010) and Zhu et al. (2012) however find that 
the relevant velocity with which they can reproduce the numer- 
ically obtained accretion rates is not the free fall velocity, but 
the Keplerian shear velocity across the Hill sphere. The veloc- 
ity with which the newly accreted matter joins the clump is thus 
of order £1Rh where Q is the Keplerian frequency. In the exam- 
ple, this corresponds to about 0.6 km/s. Whether or not this is 
sufficient to produce a supercritical shock must be evaluated for 
realistic velocities, pre-shock densities and temperatures. If it is 
not, as we may suspect due to the rather low velocity, then the 
material accreted in this phase should be added at a high spe- 
cific entropy. This leads to the second supposition: If a GI planet 
accretes most of its final mass before it undergoes the second 
collapse, and the accretion onto the AU sized clump is relatively 
gentle, we may again expect a "hot start". 

In the second sub-scenario, the clump accretes most of its 
final mass after it has undergone the second collapse which 
leads to an object with a size of a few Jovian radii (Helled & 
Bodenheimer 2011). Now, the accreted gas hits the protoplanet 
at a much higher velocity of several 10 km/s as in the runaway 
gas accretion phase in the core accretion scenario. This leads to 
the third supposition: The larger the fraction of the final mass in 
a GI planet accreted after the clump has undergone the second 
collapse, the "colder" it should be, provided that the accretion 



shock in the post-collapse phase is supercritical. The effect that 
the higher the fraction of the final mass processed through the 
shock, the lower the initial luminosity, is seen in the simulations 
of Marley et al. (2007) where it causes more massive planets 
to be less luminous than their lighter counterparts ("luminosity 
inversion"). 

This discussion (see also the related discussion in Spiegel & 
Burrows 2012) shows that with the currently limited knowledge 
both about the specific way matter is accreted, and the involved 
thermodynamics in general, it is difficult to firmly associate a 
"hot/cold start" with a specific formation mechanism. Some of 
the general trends for the impact of formation on the initial lumi- 
nosity of core accretion planets will be presented in a dedicated 
paper (Mordasini et al. in prep.). 

8.2. "Hot start" and "cold start" models 

In the present paper, we will only compare our results with al- 
ready published studies, similar as for the radii in the previous 
section. Figure 10 shows the (total) luminosity (including the ac- 
cretional luminosity) as a function of time for planets of 1, 2, 5 
and 10 M%. 

The left panel shows the result using "cold start" boundary 
conditions, i.e. it is assumed that all accretional energy is ra- 
diated away at the shock, while no losses are assumed to oc- 
cur in the "hot start (accreting)" case (Eq. 13). The cold start, 
M = \M\ simulation is the same one as discussed in detail in 
Sect. 6, while the other masses correspond to the simulations 
shown in Fig. 6. We plot in the figures also the L(t) found by 
Burrows et al. (1997) and Baraffe et al. (2003), for comparison 4 . 
Both these models are classical "hot start" simulations. We asso- 
ciate the "t = 0" moment of these models with the moment when 
in our simulations, the planets reach their final mass. For the 
"cold start" simulations, this moment corresponds to the sharp 
drop of L particularly well visible for the 10 and 5 M% cases, 
when L acc vanishes. 

8.2.1. "Hot start": comparison with Burrows et al. (1997) and 
Baraffe et al. (2003) 

Focussing first on the right panel with the "hot start (accreting)" 
models, we see a good agreement between our model and the 
two other ones. The differences between our (simpler, grey atmo- 
sphere) model and the two more complex models is of the same 
order as the differences mutually between Burrows et al. (1997) 
and Baraffe et al. (2003). This is in agreement with Bodenheimer 
et al. (2000) who also find very good agreement of their grey at- 
mosphere models and Burrows et al. (1997). It is however clear 
that the precise shape of L(t) (there is for example a small bump 
in our cooling curves when log L/L Q x -6.25) depends directly 
on the opacities, so that we expect that our predictions are some- 
what less accurate during the long-term evolution. 

One sees that the assumption of no radiative losses at the 
shock, but still a gradual building up of the planet (as in our sim- 
ulations) leads to very similar results as the classical "hot start" 
scenario where one starts with a fully formed planet. The physi- 
cal reason for the similarity is that in both cases, no entropy sink 
exists, and that the resulting short Kelvin-Helmholtz timescales 
make that the initial conditions are quickly "forgotten" in a rapid 
convergent evolution, as discussed at the end of this section. It 
means that core accretion with a radiatively completely ineffi- 

4 The data plotted is from http://www.astro.princeton.edu/~burrows/ 
and http://perso.ens-lyon.fr/isabelle.baraffe/. 
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Fig. 10. Luminosity as a function of time for planets with a mass of 1, 2, 5 and 10 M%, as labelled in the plot. The left panel shows 
the case of "cold start" boundary conditions where all accretional energy is radiated away at the shock, while the right panels shows 
"hot start (accreting)" models where no radiative losses occur. The dotted lines show for comparison the results of Burrows et al. 
(1997), while the dashed lines show BarafFe et al. (2003). Both these models use the classical "hot start" scenario. The "cold start" 
1 M% simulation is the same as shown in Fig. 3. 



cient shock leads to very similar planetary luminosities at young 
ages as in classical "hot start" simulations. 

This result certainly underlines the importance of a detailed 
description of the shock structure, as already pointed out by 
Marley et al. (2007), because we see that the shock structure 
has an influence which is as important as the formation mech- 
anism. Progress on the shock structure is an important task for 
future studies (see Commercont et al. 201 1 for a recent study of 
the stellar case). 

Note that the assumption of dl/dr = in the envelope for 
the "hot start" scenario could lead to significant departures from 
the actual interior structure and the associated luminosity of the 
planets during the gas runaway accretion phase. This is due to 
the fact that we cannot catch the effects of a possible deep ra- 
diative zone, as discussed in Sect. 3.6. This would also affect the 
luminosities at the beginning of the evolutionary phase i.e. just 
after the end of formation. 

The comparison of the luminosity at this moment as found 
in our simulations, and as assumed in Burrows et al. (1997) and 
Baraffe et al. (2003) indicates that we have luminosities which 
are similar (for the 10M^ case) or higher (for lower masses). 
A higher luminosity means that the Kelvin-Helmholtz timescale 
is very short, so that the models converge on a timescale of only 
10 5 - 10 6 yrs. Such a rapid convergence of "hot start" and "hotter 
start" models has been found by BarafFe et al. (2002, 2009) (see 
also Marley et al. 2007 for a discussion). The effect of a variable 
luminosity on "hot start" objects, calculated as described in Sect. 
3.4.1, will be studied in future work. This should allow to better 
describe the properties of very young "hot start" objects. 



8.2.2. "Cold start": comparison with Marley et al. (2007) 

Comparing the left and the right panel makes it immediately 
clear that we recover the important findings of Marley et 
al. (2007) (which were already visible in the simulations of 
Bodenheimer et al. 2000), that the luminosity of young planets 
formed according to the "cold start" assumption can be substan- 
tially lower, at least for high-mass planets. 

Focussing on the left panel with the "cold start" simula- 
tions we see that we recover the most important effects found 
by Marley et al. (2007): we also see a "luminosity inversion" 
which means that the post-accretional luminosity (i.e. the lu- 
minosity shortly after the final mass of the planet is reached) 
is the highest for the lowest mass planet. The picture is some- 
what complicated by the fact that the low-mass planets start with 
the highest L but also cool the quickest, as their KH-timescale 
GM 2 /{RL) is the shortest. The reason for the "luminosity inver- 
sion" is that the higher the total mass, the higher the fraction of 
gas in the envelope that has been process through the entropy- 
reducing shock. The total mass of the planets at the moment 
when they start to collapse is about 87 M®, with a core mass of 
22 M @ and an envelope of 65 M e (cf. Table A. 2). This is equal 
for all planets (their evolution is identical up to the point where 
for the lower mass planets, accretion is shut off). This means 
that for a IM% planet, about 20% of its final envelope mass is 
accreted in the attached phase, and thus without going through 
the entropy-reducing shock. The ratio is proportionally smaller 
for higher mass planets, leading to the inversion. 

Another identical result is that the higher the mass, the longer 
it takes until "hot start" and "cold start" luminosities become 
equal, a result that can be understood from the KH-timescales. 
For the 10 M% planet, it takes several 100 Myrs, similar to 
Marley et al. (2007). 
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Quantitatively, there are however also some differences: we 
note that we find post-accretional luminosities for the massive (5 
and 10 M\) planets that are of order log L/L Q w -4.1 to -4.3, 
while Marley et al. (2007) rather find log L/L e * -5.7, i.e. quite 
fainter. The reason for this difference will be investigated in de- 
tail in Mordasini et al. (in prep.). 

9. Summary 

We have extended our formation model to a self-consistently 
coupled formation and evolution model for planets with a pri- 
mordial H2/He envelope. In this paper we have described the 
modification we made to the computational module that calcu- 
lates the internal structure in order to make this extension pos- 
sible. We then compared the results concerning planetary evo- 
lution with existing work. We found good agreement with more 
complex models. In a companion paper we describe further ex- 
tensions and improvements relevant during the evolution of the 
planets, like an internal structure model for the solid core as- 
suming a differentiated planet which gives realistic radii also for 
planets without significant atmospheres, or the radiogenic lumi- 
nosity of the planet's mantle. In the companion paper we also 
discuss extensively our results on planetary radii as obtained in 
population synthesis calculations. 

9.1. Jupiter: coupled in situ formation and evolution 

We have simulated the combined formation and evolution of 
Jupiter in the framework of classical core accretion models 
without migration and disk evolution (Pollack et al. 1996). We 
have shown that the upgraded model with the new simple and 
rapid method to calculate evolutionary sequences leads to results 
which are in very good agreement with those of several other 
models (Lissauer et al. 2009; Burrows et al. 1997; Baraffe et al. 
2003). We have found that the nominal model for Jupiter leads 
to the formation of a planet which has properties at 4.6 Gyrs 
which are in excellent agreement with observed values in terms 
of internal composition, radius, luminosity and surface pressure- 
temperature profile. 

9.2. Planetary radii 

We have studied the mass-radius diagram of giant planets and 
the influence of the core mass on the total radius. We have found 
that the upgraded models yields planetary radii of giant plan- 
ets which are in very good agreement with more sophisticated 
models (Fortney et al. 2007; Baraffe et al. 2008). We must as- 
sume that our cooling curves are probably less accurate for low- 
mass, core-dominated planets, possibly with heavily enriched 
envelopes. Still we show in the companion Paper II that the mod- 
els yields also for a completely different type of planet, namely a 
close-in, 4.2M ffl super-Earth planet with a tenuous ~ 1 % H.2/He 
envelope a radius which agrees with the more detailed simula- 
tions of Rogers et al. (2011) to better than 10%, corresponding 
to about 0.2 R e . We thus can calculate planetary radii for a very 
wide range of planets, which is our goal for the planet popu- 
lation synthesis calculations for which the model is eventually 
intended (see Paper II). 

9.3. Planetary luminosities 

We have used our upgraded model to study the luminosity of gi- 
ant planets of different masses as a function of time. We have 



made calculations both for the "cold start" and the "hot start (ac- 
creting)" scenario. For the "hot start (accreting)" models we have 
found cooling curves which are very similar as in Burrows et al. 
(1997) and Baraffe et al. (2003) despite the fact that we use sim- 
ple gray boundary conditions. In the "cold start" calculations we 
have recovered the result from Marley et al. (2007) that the lumi- 
nosities of massive young planets are much smaller, at least for 
the chosen parameters. In a dedicated work (Mordasini et al. in 
prep.) we revisit the luminosity of young Jupiters using a large 
suite of combined formation and evolution calculations. 

10. Conclusion 

In earlier versions of the model, we calculated the internal struc- 
ture of the gaseous envelope only during the attached, pre-gas 
runaway accretion phase. This is sufficient if one is only inter- 
ested in the final mass of the planets, and thus sufficient for com- 
parisons with planets found by the radial velocity method. Now 
we calculate the structure also during the gas runaway accre- 
tion phase (which is also the phase when the planet's radius col- 
lapses) and during the evolutionary phase at constant mass over 
Gigayear timescales. With that we now know all major quan- 
tities (mass, semimajor axis, radius, luminosity, composition) 
characterizing the planets during their entire formation and evo- 
lution. This allows to compare our population synthesis models 
directly and consistently with results coming from all major ob- 
servational techniques used to detect and characterize extrasolar 
planets (see Paper II for a comparison of the synthetic and actual 
mass-radius relationship, the predicted distribution of radii, and 
the comparison with Kepler). 

Extrasolar planet research has entered an era in which mas- 
sive amounts of observational data regarding very different sub- 
populations of planets are brought to us from different tech- 
niques. They should all be explained consistently by planet for- 
mation and evolution theory, but it is a difficult task to unite the 
different elements and observational constraints into one consis- 
tent global picture. 

We think that a fruitful approach in this situation is to work 
with theoretical models able to make testable predictions in a 
consistent way for all important observational techniques. With 
the work presented in this paper we make a development in this 
direction, allowing us to test and improve theoretical formation 
models using the wealth of data coming from observations. 
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Appendix A: Comparison calculation with Pollack 
et al. (1996) 

In this appendix, we present a detailed comparison calculation with the classical 
simulation of Pollack et al. (1996). In particular, we simulate the cases Jl, Jla, 
Jib and Jlc. Simulation Jl is the nominal case, while for the other cases, solid 
accretion gets (arbitrarily) shut off at specified times (see Figure A. 1). The con- 
sequence of this is a reduced core luminosity, which in turn leads to a faster gas 
accretion due to the reduced thermal support of the envelope, an effect described 
by Pollack et al. (1996). 

This effect is relevant in our population synthesis calculations in two situa- 
tions: during planet migration, and/or the concurrent formation of several plan- 
ets. The reason for the first situation is that with the updated (non-isothermal 



21 



C. Mordasini et al.: Characterization 

Table A.l. Settings for the Pollack et al. (1996) comparison cal- 
culations. 



case 


Jl comparison 


a [AU] 


5.2 


E s ,o [g/cm 2 ] 


10 




0.01 


1 neb L^J 


150 


^neb [dyn/cm 2 ] 


0.275 


Dust-to-gas ratio 


1/70 


Initial embryo mass [M e ] 


0.1 


Migration 


not included 


Disk evolution 


not included 


Planetesimal ejection 


not included 


Core density 


constant, 3.2 g/cm 3 


Simulation duration 


10'° yrs 


Grain opacity red. factor / opil 


1 



type I) migration model, protoplanets are often seen to first migrate outwards, 
and then back in. On their way back in, they move into the region of the planetes- 
imal disk which they have cleared from planetesimals while migrating outwards 
(see Mordasini et al. 201 la). In the second situation, one planet can migrate into 
a region which has previously been cleared by another growing planet. In both 
situations, a faster gas accretion will result. From this we see that migration and 
accretion can interact. 

With this appendix, we want to check whether our simplified luminosity 
calculation reproduces the effect found by P96, and in general compare quanti- 
tatively our results with the published ones. 

The Jl initial conditions in particular mean that the initial planetesimal sur- 
face density S p is 10 g/cm 2 as in Section 6. The grain opacities are now how- 
ever at 100 % of the interstellar value (/ opa = 1). The density of the core is 
constant at 3.2 g/cm 3 as in Pollack et al. (1996), and also the nebular bound- 
ary conditions are the same as in this work. For the J la, Jib and Jlc cases, we 
switch off planetesimal accretion at the same moments as in Pollack et al. (1996). 
The parameters used for the simulations are given in Table A. 1 . It is clear that 
the two simulation still cannot yield exactly identical results, as they differ in 
some other aspects, like a different equation of state, or a different model for the 
planetesimal-envelope interaction. 

In Figure A. 1 we show the core mass, envelope mass, and total mass for the 
four models, and in Table A. 2 we give the values of the most important quantities 
at specific moments in time. We first focus on the nominal Jl case. The crossover 
point, i.e. the moment when the core and the envelope mass are identical and 
equal to M cr is found to occur at f cr = 7.49 Myrs at M cr = 16.62M e . This 
is very close to Pollack et al. (1996) with f cr =7.58 Myrs and M cr =16.17M e . 
Note however that the very good agreement in f cr is partially a chance result, 
because it is well known that changing for example just the opacity tables or the 
EOS can lead to variations of several million years (Hubickyj et al. 2005). These 
authors found for identical conditions, but an updated version of the model used 
in Pollack et al. (1996), a f cr =6.07 Myrs and A/ cl = 16.16M ffi . At crossover, the 
core and gas accretion rate are 2.37x 10~ 6 and 1.19x 10~ 5 M ffi /yr, which are both 
in very good agreement with Pollack et al. (1996) . 

The limiting gas accretion rate of 0.01 M e /yr is reached at 8.02 Myrs, about 
half a million year after crossover. The core mass is then about 22.90 M ffi and the 
envelope mass 75.8 M s . In Pollack et al. (1996), a similar, but somewhat smaller 
gas accretion rate is reached 0.42 Myrs after crossover, at Mz = 21.5M ffi and 
Mxy = 64AM S . In contrast to Pollack et al. (1996), we have taken the calcu- 
lation also through the detached end evolutionary phase. The second luminosity 
peak occurs shortly after the limiting gas accretion is reached, at 8.036 Myrs, 
when the planet has already contracted to a radius of 4.49 R q^, and grown to 
a mass of 267 M s . The luminosity is then log(L/L Q ) = -2.77. Lissauer et al. 
(2009) find peak luminosities for planets accreting with the same limiting gas 
accretion rate of log(L/L G ) ~ -2.4 to -2.3. 

The case Jla (red lines) is identical to case Jl, except that at t = 1.5 Myrs, 
the solid accretion rate is artificially ramped down to zero on a short timescale as 
in Pollack et al. (1996). This results in a faster growth of the envelope, i.e. an in- 
crease of the gas accretion rate. Our model thus reproduces this important effect. 
In this simulation, crossover occurs already at f cr =3.38 Myrs and M cr =12.98M e . 
In Pollack et al. (1996), the values are 3.32 Myrs, and 12.24 M e which shows 
that also quantitatively, there is very good agreement. The results for the other 
two cases, Jib and Jlc, are comparable, even though the differences in timescales 
are somewhat larger. 
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Fig. A.l. Total mass (solid line), envelope mass (dotted line) and 
core mass (dashed line) for the cases Jl (black), Jla (red), Jib 
(blue) and Jlc (green). The X-symbols indicate the moment in 
time when Mz is artificially ramped down for the latter three 
simulations. 



In summary we see that there is also quantitatively a very good agreement 
between the Pollack et al. (1996) simulations and the ones presented here. 
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Table A.2. Models for the in situ formation of Jupiter. Jl and 
Jla are the comparison calculations with Pollack et al. (1996). 
The third column is our own nominal model for Jupiter's in situ 
formation presented in Section 6. It should be more realistic be- 
cause it includes a reduced grain opacity, a variable core density 
and the ejection of planetesimals. 



Phase 




Jl 


Jla 


nominal 


First luminosity peak 


Time 


0.242 


0.242 


0.188 




M 


8.13 


8.13 


8.20 




Mz 


8.12 


8.12 


8.17 




M X Y 


0.01 


0.01 


0.03 




Mz 


1.00xl0~" 


1.00x10-" 


1.40x10-" 




R 


25.8 


25.8 


26.9 




log L/L 


-4.99 


-4.99 


-4.81 


Cross over point 


Time 


7.49 


3.38 


0.817 




M 


33.23 


25.96 


32.53 




M„ 


16.62 


12.98 


16.27 




M z 


2.37X10- 6 





UOxHr 5 




Mxy 


1.19xl0~ 5 


2.14xl0~ 5 


8.33xl0 -5 




R 


61.2 


53.5 


60.4 




log L/Lq 


-6.31 


-6.84 


-5.46 


Onset of limited 


Time 


8.02 


4.125 


0.923 


gas accretion 


M 


98.7 


107.5 


86.7 




M z 


22.9 


13.0 


21.8 




Mxy 


75.8 


94.5 


64.9 




M z 


5.94xl0~" 





6.13x10-" 




R 


104.0 


108.2 


98.1 




log L/L 


-4.03 


-4.10 


-4.03 


Second luminosity peak 


Time 


8.036 


4.141 


0.940 




M 


266.9 


266.0 


264.4 




M z 


32.7 


13.0 


30.9 




M XY 


234.3 


252.9 


233.5 




M z 


4.84x10-" 





3.60x10-" 




R 


4.49 


2.32 


3.77 




log L/Lq 


-2.77 


-2.54 


-2.65 






5.25X10 4 


5.01x10" 


6.08x10" 


Present day (4.6 Gyr) 


M 


318.7 


317.8 


316.6 




Mz 


35.4 


13.0 


32.9 




Mxy 


283.3 


304.8 


283.7 




R 


1.14 


1.15 


0.99 




Pcore 


3.2 


3.2 


14.31 




log L/L 


-8.83 


-8.79 


-9.01 




L/L % 


1.73 


1.88 


1.13 




T 


128.4 


129.2 


126.3 




7"ccnl 


2.16X10 4 


2.28x1 4 


1.76x10" 



The following units are used: Time in Myrs; Masses in M e ; Accretion rates in M e /yr; Radii 
in Jovian equatorial radii equal 7.15 X 10 9 cm; Core densities in g/cm 3 ; Temperatures in K; 
Lr^ is the present day Jovian intrinsic luminosity = 8.7 X 10~ 10 Lo. 



D'Angelo, G. & Lubow, S. H. 2008, ApJ, 685, 560 

D'Angelo, G., Durisen, R. H., & Lissauer, J. J. 2010, in Exoplanets, Space 
Science Series, ed. S. Seager (University of Arizona Press, Tucson), 319 

Dodson-Robinson, S. E., Veras, D., Ford, E. B., & Beichman, C. A. 2009, ApJ, 
707, 79 

Forgan, D., & Rice, K. 201 1, MNRAS, 417, 1928 

Fortier, A., Benvenuto, O. G., and Brunini, A. 2007, A&A, 473, 311 

Fortney, J. J., Marley, M. S., Hubickyj, O., Bodenheimer, R H., & Lissauer, J. J. 

2005, Astron. Nachrichten, 326, 925. 
Fortney, J. J., Marley, M. S., & Barnes, J. W. 2007, ApJ, 659, 1661 
Fortney, J. J., & Nettelmann, N. 2010, Space Sci.Rev., 152, 423 
Fortney, J. J., Baraft'e, I., & Militzer, B. 2010, in Exoplanets, Space Science 

Series, ed. S. Seager (University of Arizona Press, Tucson), 397 
Fouchet, L., Alibert, Y., Mordasini, C, & Benz, W. 201 1, A&A, 540, A107 
Freedman, R. S., Marley, M. S., & Lodders, K. 2008, ApJS, 174, 504 
Galvagni, M., Hayfield, T., Boley, A., Mayer, L., Roskar, R. & Saha, P. 2012, 

MNRAS, submitted 
Gillon, M., Doyle, A. P., Lendl, M., et al. 2011, A&A, 533, A88 
Guillot, T. 1999, Science, 286, 72 
Guillot, T. 2005, ARA&A, 33, 493 

Guillot, T., Santos, N. C., Pont, F, et al. 2006, A&A, 453, L21 

Guillot, T. & Gautier, D. 2009, in Treatise on Geophysics, ed. G. Shubert & T. 

Spohn, 10, 439 
Guillot, T. 2010, A&A, 520, 27 

Hartmann, L., Cassen, P., & Kenyon, S. J. 1997, ApJ, 475, 770 

Helled, R., & Bodenheimer, P. 201 1, Icarus, 21 1, 939 

Heng, K., Hayek, W., Pont, R, & Sing, D. K. 2012, MNRAS, 420, 20 



Hubbard, W. B., Guillot, T, Marley, M. S., et al. 1999, Planet. Space Sci., 47, 
1175 

Hubickyj, O., Bodenheimer, P. H., & Lissauer, J. J. 2005, Icarus, 179, 415 

Ida, S. & Lin, D.N.C. 2004, ApJ, 604, 388 

Janson, M., Bonavita, M., Klahr, H., et al. 201 1, ApJ, 736, 89 

Kasper, M. E., Beuzit, J.-L., Verinaud, C, et al. 2008, Proc. SPIE, 7015 

Kippenhahn & Weigert 1990, Stellar Structure and Evolution (Springer, Berlin) 

Kratter, K. M., Murray-Clay, R. A., & Youdin, A. N. 2010, ApJ, 710, 1375 

Leconte, J., Chabrier, G., Baraffe, I., & Levrard, B. 201 1, Proceedings Detection 

and Dynamics of Transiting Exoplanets, ed. F. Bouchy, R. Diaz & 

C.Moutou, 11,03004 
Leconte, J., & Chabrier, G. 2012, A&A, 540, A20 
Leger, A., Rouan, D., Schneider, J., et al. 2009, A&A, 506, 287 
Lissauer, J. J., Hubickyj, O., D'Angelo, G, & Bodenheimer, P. H. 2009, Icarus, 

199, 338 

Lodders, K., & Fegley, B. 1998, The planetary scientist's companion (Oxford 

University Press, Oxford) 
Lubow, S. H, Seibert, M., & Artymowicz, P. 1999, ApJ, 526, 1001 
Lubow, S. H. & D'Angelo, G. 2006, ApJ, 641, 526 
Marois, C, Macintosh, B., Barman, T, et al. 2008, Science, 322, 1348 
Marois, C, Zuckerman, B., Konopacky, Q. M., Macintosh, B., & Barman, T. 

2010, Nature, 468, 1080 
Marley, M. S., Fortney, J. J., Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 

2007, ApJ, 655, 541 

Marley, M. S., Saumon, D., Cushing, M., et al. 2012, ApJ, in press 
(arXiv: 1205.6488) 

Mayor, M., Marmier, M., Lovis, C, et al. 2011, A&A, submitted 

(arXiv: 1109.2497) 
McBride, J.. Graham, J. R., Macintosh, B., et al. 201 1, PASP, 123, 692 
Mercer-Smith, J. A., Cameron, A. G. W., & Epstein, R. I. 1984, ApJ, 279, 363 
Militzer, B., Hubbard, W. B., Vorberger, J., Tamblyn, I., & Bonev, S. A. 2008, 

ApJ, 688, L45 
Miller, N., & Fortney, J. J. 201 1, ApJ, 736, L29 

Mizuno, H., Nakazawa, K., & Hayashi, C. 1978, Prog. Th. Phys., 60, 3, 699 
Mordasini, C, Alibert, Y, & Benz, W. 2005, in Proceedings of Haute Provence 

Observatory Colloquium, ed. L. Arnold, 85 
Mordasini, C, Alibert, Y, & Benz, W. 2009a, A&A, 501, 1 139 
Mordasini, C, Alibert, Y, Benz, W., & Naef, D. 2009b, A&A, 501, 1 161 
Mordasini, C, Dittkrist, K. M., Alibert, Y, et al. 2011a, IAU Symposium, 276, 

72 

Mordasini, C, Mayor, M., Udry, S., Lovis, C, Segransan, D., Benz, W., et al. 

201 lb, A&A, 526, 111 
Mordasini, C, Alibert, Y, Klahr, H., & Benz, W. 2011c, Proceedings Detection 

and Dynamics of Transiting Exoplanets, ed. F. Bouchy, R. Diaz & 

C.Moutou, 11,04001 
Mordasini, C, Alibert, Y, Benz, W., Klahr, H., & Henning, T. 2012a, A&A, 

541, A97 

Mordasini, Alibert, Y, Georgy, C, Dittkrist, K.-M., Henning, T. 2012b, A&A, 

in review (Paper II) (arXiv: 1206.3303) 
Mordasini, C, Klahr, H., Alibert, Y, Henning, T. & Benz, W. 2012c, A&A, in 

review 

Movshovitz, N., Bodenheimer, P. H., Podolak, M., & Lissauer, J. J. 2010, Icarus, 
209,616 

Nettelmann, N., Becker, A., Hoist, B., & Redmer, R. 2012, ApJ, 750, 52 

Papaloizou, J. C. B. & Terquem, C. 1999, ApJ, 521, 823 

Papaloizou, J. C. B. & Nelson, R. P. 2005, A&A, 433, 247 

Perri, F, & Cameron, A. G. W. 1974, Icarus, 22, 416 

Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62 

Prialnik, D., & Livio, M. 1985, MNRAS, 216, 37 

Rafikov, R. R. 2007, ApJ, 662, 642 

Rogers, L. A., Bodenheimer, P., Lissauer, J. J., & Seager, S. 201 1, ApJ, 738, 59 

Saumon, D., Chabrier, G. & Van Horn, H. M. 1995, ApJS, 99, 713 

Saumon, D., & Guillot, T. 2004, ApJ, 609, 1 170 

Shiraishi, M. & Ida, S. 2008, ApJ, 684, 1416 

Showman, A. P., & Guillot, T. 2002, A&A, 385, 166 

Spiegel, D. S., & Burrows, A. 2012, ApJ, 745, 174 

Stahler, S. W., Shu, R, & Taam, R. E. 1980, ApJ, 241, 637 

Stevenson, D. J. 1979, MNRAS, 187, 129 

Stevenson, D. J. 1985, Icarus, 62, 4 

Valencia, D., Ikoma, M., Guillot, T, & Nettelmann, N. 2010, A&A, 516, A20 

Wilson, H. R, & Militzer, B. 2012, Physical Review Letters, 108, 111101 

Wolfgang, A., & Laughlin, G. 2012, ApJ, 750, 148 

Wuchterl, G. 1991, Icarus, 91, 53 

Zapolsky, H. S. & Salpeter, E. E. 1969, ApJ, 158, 809 

Zhou, J.-L. & Lin, D. N. C. 2007, ApJ, 666, 447 

Zhu, Z., Hartmann, L., Nelson, R. P., & Gammie, C. F. 2012, ApJ, 746, 110 



23 



